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' ABSTRACT 

We investigate the importance of several numerical artifacts such as lack of resolu- 
tion on spectral properties of the Lyman-a -forest as computed from cosmological 
hydrodynamic simulations in a standard cold dark matter universe. We use a new 
simulation code which is based on a combination of a hierarchical particle-particle- 
particle-mesh (P3M) scheme for gravity and smoothed particle hydrodynamics (SPH) 
^ ' for gas dynamics. We have performed extensive comparisons between this new code 

Q> I and a modified version of the HYDRA code of Couchman et al. and find excellent 

agreement. We have also rerun the TREESPH simulations of Hernquist et al. using 
our new codes and find very good agreement with their published results. This shows 
\ that results from hydro dynamical simulations that include cooling are reproducible 

with different numerical algorithms. We then use our new code to investigate several 
numerical effects such as resolution on spectral statistics deduced from Voigt profile 
fitting of lines by running simulations with gas particle masses of 1.4 x 10^, 1.8 x 10^, 

■ 2.2 X 10^ and 2.1 x IQ^Mq. When we increase the numerical resolution the mean ef- 
PU| fective hydrogen optical depth converges and so does the column density distribution. 

However, higher resolution simulations produce narrower lines and consequently the 
• 6-parameter (velocity width) distribution has only marginally converged in our high- 

^ \ est resolution run. Obtaining numerical convergence for the mean Heii transmission 

. is demanding. When progressively smaller halos are resolved at better resolution, a 

larger fraction of low density gas contracts to moderate over densities in which Hen is 
already optically thick, and this increases the net transmission, making it difficult to 

■ simulate Hen reliably. Our highest resolution simulation gives a mean effective optical 
5-H ' depth in Hen 5% lower than the simulation with eight times lower mass resolution, 

illustrating the degree to which the Hen optical depth has converged. In contrast, the 
hydrogen mean optical depth for these runs is identical. Since many properties of the 
simulated Lyman-a -forest depend on resolution, one should be careful when deduc- 
ing physical parameters from a comparison of the simulated forest with the observed 
one. We compare predictions from our highest resolution simulation in a cold dark 
matter universe with a photo-ionising background inferred from quasars as computed 
by Haardt & Madau (1996), with observations. The simulation reproduces both the 
Hi column density and 6-parameter distribution when we assume a high baryon den- 
sity, flsh^ ~ 0.028. In addition we need to impose a higher IGM temperature than 
predicted within our basic set of assumptions. We argue that such a higher temperature 
could be due to differences between the assumed and true reionization history. The 
simulated Hi optical depth is in good agreement with observations but the Hen optical 
depth is lower than observed. Fitting the Hen optical depth requires a larger jump 
~ 14 between the photon flux at the Hi and Hen edge than is present in the Haardt 
& Madau spectrum. 

Key words: cosmology: theory - hydrodynamics - large-scale structure of universe 
- quasars: absorption lines 
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1 INTRODUCTION 

Sight lines to distant quasars intersect many cosmological 
structures containing neutral hydrogen and Lyman-a scat- 
tering by the Hi in these structures produces a forest of 
lines blueward of the quasar's Lyman-a emission line (Lynds 
1971). This 'Lyman-Q forest' contains unbiased information 
on the temperature, density, velocity and ionization struc- 
ture of the intergalactic medium (IGM) along the line of 
sight to the quasar, making the structures responsible for 
the Lyman-a forest a useful probe for studying the high- 
redshift universe. In addition, it is likely that the absorbing 
gas retains a memory of its state at even higher redshifts, 
enabling us to study its initial conditions (Croft et al. 1998) 
and previous history. Since these structures are of moder- 
ate density contrast, they are easier to simulate numerically 
than galaxies, and consequently the high redshift universe 
can be studied efficiently and accurately by comparing sim- 
ulations of the Lyman-a forest with observations. 

Recent hydrodynamic simulations of hierarchical struc- 
ture formation in a universe dominated by cold dark matter 
(CDM) have been shown to be remarkably successful in re- 
producing a variety of statistics of Lyman-a absorption lines 
(Cen et al. 1994, Zhang, Anninos & Norman 1995, Miralda- 
Escude et al. 1996, Hernquist et al. 1996, Wadsley & Bond 
1996, Zhang et al. 1997), including the number of lines per 
unit redshift per unit column density and the number of lines 
with given width ('6' parameter), as well as its evolution 
at low-redshift (Theuns, Leonard & Efstathiou 1998). This 
is quite encouraging for the hierarchical picture of structure 
formation since the underlying cosmological models were de- 
signed with galaxy formation in mind, hence their Lyman- 
a properties can be considered to be a genuine and successful 
prediction. Most simulations to date have assumed a critical 
density, cold dark matter model, in which a photo-ionising 
background close to that inferred from quasars as computed 
by Haardt & Madau (1996) is required to explain the prop- 
erties of the Lyman-Q forest. However, other variants of the 
CDM model still provide acceptable fits, with only minor 
modifications to the required photo-ionization background 
(Cen et al. 1994, Miralda-Escude et al. 1996). 

In this paper we introduce a new simulation code de- 
signed to study numerically the formation of Lyman-o? sys- 
tems. It is based on a combination of Smoothed Particle Hy- 
drodynamics (SPH, Lucy 1977, Gingold & Monaghan 1977, 
see e.g. Monaghan 1992 for a review) and an adaptive P3M 
(particle-particle-particle-mesh) gravity solver (Couchman 
1991). Its efficient gravity solver and SPH implementation 
lead to a fast and accurate code which has the potential to 
extend considerably the dynamic range of the simulations. 
We discuss tests of the new code and perform extensive 
comparisons against two other simulation codes: HYDRA 
and TREESPH . Both of these are also based on SPH but 



their g ravity solvers differ: HYDRA (Couchman, Thomas & 
Pearce 1995) uses the same gravity solver as APMSPH but 



TREESPH ( |Katz, Weinberg fc Hernquist 1996b| ) uses a tree 
structure. We discuss in detail the differences between the 
APMSPH and HYDRA codes. We also discuss the changes 
we have made to the publically available HYDRA code to 
study the Lyman-Q cloud problem. The overall agreement 
between the three codes is excellent which shows that hydro- 
dynamic simulations that include cooling are reproducible 



with different simulation codes. The good agreement also 
shows that HYDRA can be used to study the Lyman-a prob- 
lem and we are currently analysing several large HYDRA 
simulations performed on the T3D computer to understand 
in more detail how resolution affects Lyman-Q statistics (the 
VIRGO consortium, in preparation). 

We then use APMSPH to perform simulations at in- 
creased resolution and establish the extent to which pub- 
lished results are influenced by lack of numerical resolu- 
tion and other numerical artifacts. Wadsley & Bond (1996; 
see also Bond & Wadsley 1997) recently warned simulators 
of the importance of long-wavelength perturbations on the 
occurrence of filamentary structures in simulations. This 
is illustrated explicitly in the work of Miralda-Escude et 
al. (1996) who compare simulations with the same resolu- 
tion but different box sizes. Unfortunately, current numerical 
codes do not possess the required dynamic range to resolve 
the Jeans length in a very large simulation box. We try to 
gauge the effects of missing waves and of failing to resolve 
the Jeans length by performing simulations with various box 
sizes. This paper is organised as follows: Section ^ discusses 
the physical model and gives details of the simulation codes. 
Section ^ presents the comparisons between codes. Section 
addresses the importance of numerical resolution. Section 
does a comparison of simulations against observations and 
finally Section ^ summarises. Technical details are relegated 
to Appendices. 



2 SIMULATION 
2.1 Physical model 

We model the evolution of a periodic, cubical region of a crit- 
ical density Einstein-de Sitter universe (f2 = 1, JIa = 0). The 
Newtonian equations of motion governing the evolution of 
structures are given in Appendix A. The comoving size of the 
simulated box is L/{2h) Mpc, where the Hubble constant to- 
day is written as Ho = lOO/i km s~^ Mpc~^. We will assume 
h — 0.5 throughout and describe simulations with L = 2.5, 
L = 5.5, 11.11 and 22.22Mpc. A fraction Qb = 0.05 of the 
matter density is assumed to be baryonic, consistent with 
limits from nucleosynthesis (Walker et al. 1991, but note 
the continuing debate on the deuterium abundance derived 
from quasar spectra, favouring higher values fis ~ 0.075, 
see e.g. Buries & Tytler 1997 and references therein). The 
rest of the matter is in the form of cold dark matter (DM). 
These model parameters were chosen to enable comparison 
with the TREESPH simulation of Hernquist et al. (1996) 
which has identical parameters to our lowest resolution run. 
We use the smaller boxes which have correspondingly higher 
resolution to test for numerical convergence. The simulations 
are started at a redshift z — 49 and we follow the evolution 
to z — 2. To generate initial conditions for the simulations 
we use the following fit to the CDM linear transfer function 
from Bardeen et al. (1986) 

T{k) = (l + 3.89g + (16.1g)V(5.46g)^-f (6.719)*)-'/" 
ln(l + 2.34g) 



2.34g 



(1) 



where q = k/h^Mpc and normalise it such that as = 0.7. 
Here, (t| denotes the r.m.s. of mass fluctuation in spheres 
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of radius 8h ^ Mpc. This value of ag is higher than the 
one deduced from the abundance of galaxy clusters (ag = 
0.52 ± 0.04, Eke, Cole & Frenk 1996) but we use it to allow 
a direct comparison with the Hernquist et al. (1996) simu- 
lations. The equations of motion used in the SPH code are 
detailed in Appendix Al. 



Gas is allowed to interact with the cosmic microwave 
background radiation (CMB) through Compton cooling and 
with an imposed uniform background of ionising photons, 
assumed to originate from quasars and/or young galaxies. 
Collisions between atoms that lead to ionization represent a 
loss term for the optically thin gas causing cooling, whereas 
photo-ionizations heat the gas because the ionised electron 
carries excess kinetic energy. We detail the temperature de- 
pendence of the cross-sections for these processes in Ap- 
pendix ^. They are taken from Cen (1992) with some mi- 
nor modifications. The fiux spectrum of the ionising pho- 
tons from a given quasar source seen by an average Lyman- 
Q cloud is changed due to reprocessing (absorption and emis- 
sion) by intervening clouds in the clumpy IGM. The ampli- 
tude of the ionising background changes due to the evolution 
of the quasar luminosity function, causing the flux spectrum 
to depend on redshift. Haardt & Madau (1996) took all these 
effects into account and provide fits to the photo-ionization 
and photo-heating rate as a function of redshift for the case 
where the main sources of UV photons are quasars. We use 
the Haardt & Madau rates but divide them by two so that 
they are identical to the rates used in the TREESPH sim- 
ulations of Hernquist et al. (1996), enabling us to compare 
directly with their results. The imposed background flux is 
time- dependent but we assume nevertheless that the gas re- 
mains in ionization equilibrium throughout (see below) . The 
resulting fits to the photo-heating and photo-ionization rates 
can be found in Appendix ^ and we will refer to them as 
'HM/2'. We will also present simulations with an ionising 
background of radiation with constant amplitude J21 and 
power law spectrum with index a (see Appendix ^ for de- 
tails). J21 denotes the amplitude of the radiation spectrum 
at the hydrogen Lyman-a edge in the usual units (lO"^'^ ergs 
cm~^ s~^ Hz"'^ sr~^). The net cooling and heating rates de- 
pend on the relative helium abundance by mass for which 
we assume Y = 0.24. 

Since our simulations include gas we need to specify 
the initial temperature of the IGM at high redshift. We will 
first describe the temperature evolution in the absence of 
any extra energy input and then discuss various reioniza- 
tion scenarios. These models are shown in figure |l]and were 
computed using a scheme based on Giroux & Shapiro (1996). 
At high redshifts, Compton cooling is very efficient hence 
gas will cool very quickly due to the coupling of the free 
electrons left over from recombination to the CMB photons 
giving T ~ TcMB = 2.7 x (1 -|- z)K. At later times, the gas 
becomes progressively more neutral to make the coupling 
with the CMB inefficient, causing the temperature to drop 
almost adiabatically, T oc (1 -I- z)" , with a ~ 1.8 (Giroux & 
Shapiro 1996). This temperature evolution is shown as the 
dashed line in figure 0. In contrast, in the simulations de- 
scribed here, we assume an IGM temperature at the mean 
IGM density of T = lO^K at 2 = 49 and in addition as- 
sume ionization equilibrium, in which case the gas will cool 
adiabatically T oc (1 + z)'^ (full line in the upper panel). 
If the starting temperature is higher, the gas will cool very 




Figure 1. Evolutionary tracks of IGM temperature at p/ps = 1 
(n = 1, h = 0.5, = 0.05, Y = 0.24). (a) Evolution from z = 49 
assuming the background ionising flux as computed by Haardt & 
Madau (1996) but with amplitude divided by 2, ionization equi- 
librium and setting initial temperatures, = 10'* K (solid), lO'*'^ 
and lO"^'^ K (both dotted). The evolution in the non-equilibrium 
case, with Ti = 2.74(1 -|- 49) is also plotted (dashed). Both cases 
are also shown in the lower panel for comparison (dotted and 
dashed respectively), (b) Temperature evolution including non- 
equilibrium effects resulting when a constant ( J21 = 0.5, a = 1) 
power-law UV background is 'switched on' at Zon = 9, 7.5 and 6. 
These converge with decreasing 2 to the temperature where the 
net cooling time equals the Hubble time (dot-dashed). We also 
show evolution with Zon = 9 assuming o = 3, similarly converg- 
ing on the (long dashed) temperature for which the net cooling 
time equals the Hubble time for that radiation spectrum. 



efficiently through Compton cooling until it becomes mostly 
neutral at T ~ 10'*K, after which it cools adiabatically (up- 
per dotted line). If the starting temperature is lower, it will 
cool adiabatically from the start (lower dotted curve). In 
any of the previous cases, T will be low with respect to the 
reionization value provided reionization occurs at redshifts 
z < 10 say, hence T after reionization is independent of our 
assumed starting temperature at z = 49. 

The behaviour of T after reionization depends on the 
reionization scenario and to some extent on whether non- 
equilibrium ionization effects are taken into account, as is 
illustrated in figure ^ (lower panel). If reionization occurs 
impulsively, more ionizations will take place per unit time 
than in equilibrium conditions so the non-equilibrium gas 
will become hotter (Miralda-Escude & Rees 1994, Haehnelt 
& Steinmetz 1998). Since the thermal time scales are long 
in low density gas, this difference in temperature may sur- 
vive to low redshifts. The IGM temperature is sensitive 
to the imposed photo-ionization heating rate, notably the 
Hell heating rate at reionization, which is based on an un- 
certain extrapolation to 2 ~ 6 of the quasar luminosity func- 
tion. This introduces a considerable uncertainty in the tem- 
perature of the IGM, even at redshifts as low as ~ 2. 
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The temperature of the non-equihbrium gas wiU ap- 
proach that of the equihbrium gas if the flux of photo ion- 
ising photons evolves slowly after reionization so that ion- 
ization equilibrium is re-established. From then on, photo 
ionization heating competes with Compton cooling and adi- 
abatic expansion in trying to change the gas temperature. 
The age of the universe at a given redshift then basically 
determines how long the gas has been heated, which in turn 
determines its temperature. This sequence of events is illus- 
trated in figure for a range of reionization scenarios. In the 
low density regime, the gas temperature depends on its den- 
sity in a straightforward way which is derived in Appendix ^ 
(see also Giroux & Shapiro 1996 and Hui & Gnedin 1997). 
The equilibrium temperature of the IGM, i.e. that tempera- 
ture where heating balances cooling, is generally higher than 
this {e.g. Zhang et al. 1997). In contrast, at high densities 
line cooling in the shock-heated gas becomes important and 
the gas cools to temperatures of a few times lO^K, below 
which atomic hydrogen cooling becomes inefficient. Since 
we do not include molecular hydrogen, the gas in the sim- 
ulations cannot cool radiatively below T = lO^K. In sum- 
mary (for 2 ~ 6 — > 2, after reionization): for low densities 
(< 0.1{pb)) there is a one-to-one relation between T and p 
obeyed by unshocked gas which is determined by the photo 
ionization heating rate, adiabatic cooling, and to a small ex- 
tent by Compton cooling. At high densities, p > 10^ (ps), 
collisions try to cool the shock-heated gas to T ~ 10*K. At 
intermediate densities, there is a large range in T for any 
given p, depending on the extent to which the gas has been 
shocked. 



2.2 Simulation methods 

The numerical methods described here are based on 
Smoothed Particle Hydrodynamics (SPH, Lucy 1977, Gin- 
gold & Monaghan 1977, see e.g. Monaghan 1992 for a re- 
view) for hydrodynamics and P3M (Hockney & Eastwood 
1988, Efstathiou et al. 1985) for self-gravity. The gas in the 
simulation is represented by a set of SPH particles which 
each carry the same mass but possibly a different thermal 
energy. Spline interpolation over these particles allows one 
to compute smoothed estimates for density and tempera- 
tures throughout the computational volume. Gradients may 
also be computed. The width of the spline kernel is matched 
to the local particle number density and so high density 
regions have higher numerical resolution than do low den- 
sity regions, in contrast to Eulerian schemes. P3M uses a 
combination of Fast Fourier Transforms (FFTs) and local 
direct summation to combine speed with accuracy. We com- 
pare detailed results from two codes. APMSPH was written 
specifically for this problem by one of us (TT) and is based 
on the hierarchical P3M code of Couch man (1991). HYDRA 
(Gouchman, Thomas & Pearce 1995) is a publically avail- 



able code used extensively by the VIRGO consortium. We 
have modified HYDRA to include photo-heating effects and 
to improve the simulation of low density regions (Section 



2.4 



and Appendix A3). We also c ompare with the published re- 



suhs of the TREESPH code ( Hernquist fc Katz 1989 



Katz 



Weinberg & Hernquist 1996b ) . In the rest of this section we 
will describe some technical details of the codes pertinent to 
their comparison. 



2.3 APMSPH 

This code is based on the adaptive P3M implementation of 
Couchman (1991) which uses mesh refinements in regions of 
high particle number density to speed-up gravity particle- 
particle interactions. The SPH implementation uses a sim- 
ilar but separate linked-list scheme based on a hierarchy of 
grids to find neighbouring SPH particles. This is because 
the refinements used by the gravity part of the code are not 
necessarily optimal for the SPH calculation and vice versa. 
Note that in APMSPH we find the neighbours of all particles 
even in low density regions, in contrast to other P3M-I-SPH 
implementations such as HYDRA (see below) and Evrard's 
(1988) versions. The explicit expressions used to compute 
the SPH accelerations are given in Section AS as well as 
details of our method of determining SPH neighbours. 

Given a power spectrum we set-up initial conditions 
by perturbing particles from a grid using the Zel'dovich 
(1970) approximation (Efstathiou et al. 1985). We take the 
DM and SPH particle grids off'set by half a cell size. We 
then march the particle positions forwards in time using 
a second-order accurate leap-frog integrator with variable 
time-step, using the correction procedure of Hernquist & 
Katz (1989) to keep the scheme second-order even when the 
step changes. The SPH accelerations depend on the parti- 
cles' velocities: we synchronise these by predicting velocities 
over half a time-step. As is usual, we take a time-step based 
on the Courant-Friedrichs-Levy condition (CFL, e.g. Mon- 
aghan 1992), but take smaller steps whenever violent shocks 
occur. Such shocks are fiagged by large values of the artifi- 
cial viscosity terms {e.g. Katz et al. 1996b) and the latter 
decrease the allowed time-step for accuracy and numerical 
stability {e.g. Hockney & Eastwood 1988, ch. 4). Since we 
use a uniform time-step, we take as system time-step the 
minimum time-step over all particles. To avoid one or a few 
particles from unnecessarily slowing down all the others, we 
make sure that there are a reasonable number of particles 
with similarly small time-step, by increasing the resolution 
length of those few offending particles that would otherwise 
require an even smaller step by a factor ~ 1.2. 

The integration of the cooling terms requires pro- 
hibitively small time-steps, comparable to the local cooling 
time. We solve this problem in the usual way by integrating 
the thermal energy equation at with an implicit scheme, as- 
suming in addition ionization equilibrium and fixed density. 
We evaluate the cooling and heating rates by interpolation 
from tables which we recompute every time the background 
flux changes significantly. 

This new SPH-I-AP3M implementation was tested with 
a variety of test problems: the ID shock tube, the spherically 
symmetric collapse of a gas sphere and the formation of a 
massive galaxy clusterQ Energy conservation based on the 
the Layzer-Irvine cosmic energy (equation A£ below) gives 
typically l^I /W ~ 1 — 2% for the runs discussed here. A 
simulation with 64"^ SPH and equal number of DM particles 
takes ~ 1000 steps to evolve from z = 49 to z = 2 for a box 
size of 22Mpc. Using a 128"^ grid for the gravity calculations, 
the code takes 110s per step at z = 50 (40s for gravity, 50 for 



* We are indebted to Adrian Jenkins for providing us with the 
initial conditions and final profiles for this problem as computed 
using HYDRA 
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SPH) which increases to 230s at z = 2 (125s for gravity, 82 
for SPH) on a DEC-alpha 4100 server. The total simulation 
then takes ~ 32 hours and requires ~ 60 M bytes of RAM 
memory in single precision (32 bits). 



2.4 HYDRA 

HYDRA (Couchman et al. 1995) is similar to APMSPH in 
that it combines a variable resolution SPH code with the 
adaptive P3M code of Couchman (1991). This code has been 
used extensively by the VIRGO Consortium on massively 
parallel computers such as the Cray T3D & T3E computers 
in Edinburgh and Munich, with the primary aim of study- 
ing large scale structure (Jenkins et al. 1997), cluster (Col- 
berg et al. 1997) and galaxy formation problems (Pearce et 
al. 1998), where typically very large dynamic ranges are re- 
quired. Throughout its development the code has also been 
used to explore the accuracies and inaccuracies of the SPH 
technique in modelling astrophysical scenarios, in keeping 
with the spirit of this paper. In this section we concentrate 
on reporting the code-modifications introduced to allow an 
accurate treatment of photoionised primordial gas dynamics 
on scales of ~ lOkpc. 

The first modification involved updating the thermal 
energy solver to take into account the cooling and heating 
functions in ionization equilibrium, as listed in Appendix H 
and also used in APMSPH . We then tested this new solver 
by comparing with the APMSPH one and another solver 
developed by one of us (APBL, based on a scheme of Giroux 
& Shapiro 1996), for a wide range of initial temperatures, 
densities and photo ionising background choices. 

Since the Lyman-a forest absorption is expected to arise 
from gas in the IGM at moderate over densities and slight 
under densities, it is paramount that an accurate handling 
of low-density gas is achieved in the simulations. HYDRA 
suffers in this respect since its implementation requires that 
the same linked-list to find neighbours be used for both grav- 
ity and SPH interactions. This leads to the following prob- 
lem: if a linked-list is used that is optimal for the gravity 
calculation, then particles in low density regions will not 
find enough SPH neighbours to accurately compute hydro- 
dynamic quantities such as their density. If, on the other 
hand, a linked-list is used which is suited to the SPH com- 
putation, then the gravity calculation will become very time- 
consuming. We decided to deal with this problem in a two- 
step fashion. Firstly, we run our simulations with FFT cells 
for the base mesh (the PP linking cells are typically 2-3 FFT 
cells on a side) set at half the mean inter-particle separation 
(leading to an efficient gravity calculation) but use a cor- 
rection for the SPH smoothing kernel applied for particles 
with few neighbours (see below). Secondly, post-simulation, 
we continue every end state for two further time-steps where 
the number of linking cells is decreased by a factor ^ and 
the original SPH kernel is used to re-calculate the densities 
for all particles with p/pB ~ 10 from their respective posi- 
tions as simulated. Only particle densities are ever updated 
in this step, referred to as the 'density reconstruction' step. 

Since the temperatures in the low-density IGM are den- 
sity dependent these must be post-adjusted during density 
reconstruction as well. Fortunately we can use the power- 
law density temperature relation obeyed by low density gas 



(see figure A2 below). This relation is history independent if 



reionization occurs grad uall y or at sufficiently high redshift 
as discussed in Section 



2.1 



In Section ^ we shall compare 
the effects of this scheme with the exact density scheme of 
APMSPH , and show that the temperature-density distri- 
butions simulated by the two codes are essentially indistin- 
guishable, also at low densities. More details of this HYDRA 



low density correction are given in Section A3 



2.5 TREESPH 



TREESPH (Katz, Weinberg fc Hernquist 1996b) uses a tree 



structure to compute gravitational forces using multipole ex- 
pansion based on an accuracy criteria and to find SPH neigh- 
bours. Each SPH particle is advanced with its own time-step 
based on the CFL condition, which may be smaller than the 
' system' time-step imposed on all the coUisionless parti- 
cles. TREESPH uses the same spline kernel as APMSPH 
and HYDRA but a different symmetrization procedure to 
compute SPH pairwise quantities (Katz et al. 1996b). Grav- 
itational forces are softened using a spline-kernel softening 
and made periodic using the Ewald summation method. We 
will compare our results with the published TREESPH sim- 
ulation results of Hernquist et al. (1996). 



3 CODE COMPARISONS 

The characteristics of the absorption lines in simulated spec- 
tra reflect the properties of the neutral gas in the simulation 
box. For example, stronger lines tend to be produced in re- 
gions with higher densities and correspondingly higher tem- 
peratures. For such lines, effects of peculiar velocities and 
Hubble expansion are less important and deviations from 
Voigt profiles are usually caused by surrounding matter that 
often lies in a filamentary distribution and contributes to 
the absorption. Consequently, the statistical properties of 
strong lines mostly refiect the statistics and shapes of halos 
that form in the simulation. In contrast, weaker lines usually 
originate in filaments and pancakes and residual Hubble ex- 
pansion contributes significantly in shaping the line. These 
lines then refiect the properties of the gas distribution in the 
low density IGM. We will first compare the overall gas dis- 
tribution between APMSPH and HYDRA by studying the 
amount of cold gas responsible for most of the weaker lines 
before concentrating on the halo statistics. In the second 
part of the section, we will compute and compare simulated 
spectra. These will be analysed using Voigt profile fitting in 
the same way as observed spectra. 

To make comparison with the TREESPH runs (Hern- 
quist et al. 1996) we have run simulations with initial condi- 
tions computed from the initial linear density field originally 
used by these authors which they kindly provided for us. We 
shall refer to these runs as 22-64-k, where 22 refers to the box 
size (22.22 comoving Mpc) and 64 to the cube root of the 
number of particles {i.e. 64'' particles of a single species) in 
the simulation. Throughout this paper we will compare re- 
sults from simulations using a variety of different codes and 
numerical parameters and refer to them as above. These are 
summarised in Table |l| and are all performed in a standard 
CDM cosmology {Q. = 1, h = 0.5). In what follows, APM- 
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Table 1. Runs performed. All runs assume a standard adiabatic, scale-invariant CDM cosmology = 1, 
f^A = Oi ^'^O = 50 km s~^ Mpc~^) with normalisation crs = 0.7 and baryon fraction = 0.05. Each run had 
equal numbers of SPH and dark matter particles. HKWM refers to Hernquist et al. (1996) initial conditions, 
e is the comoving softening for the spline gravity kernel, the SPH resolution length h is not allowed to fall 
below e/2. 'HM/2' refers to the radiation spectrum computed by Haardt &: Madau (1996), with amplitude 
divided by 2. J21 = 0.5, a = 1 refers to a power law spectrum of ionising photons, cen and apb refer to 
different fits to the cooling and heating functions (see Appendix . 





Name 


Box size (Mpc) 


A^SPH 


ICs 


code 


e (kpc) 


J{z) 


rates 


1 


22-64-k 


22.22 


643 


HKWM 


APMSPH ,HYDRA 


17.3,20 


HM/2 


cen 


2 


22-64 


22.22 


643 


own 


APMSPH ,HYDRA 


17.3,20 


HM/2 


apb 


3 


11-64 


11.11 


643 


own 


APMSPH ,HYDRA 


20,10 


HM/2 


apb 


4 


5-64 


5.5 


643 


own 


APMSPH 


10 


HM/2 


apb 


5 


2.5-64 


2.5 


643 


own 


APMSPH 


4.5 


HM/20 


apb 


6 


22-32 


22.22 


323 


own 


APMSPH ,HYDRA 


17,20 


HM/2 


apb 


7 


11-32 


11.11 


323 


own 


APMSPH ,HYDRA 


8.6,10 


HM/2 


apb 


8 


5-32 


5.5 


323 


own 


APMSPH 


4.3 


HM/2 


apb 


9 


11-64-j 


11.11 


643 


own 


HYDRA 


10 


J21 = 0.5, a = 1 


apb 



" For this run we imposed the 2 = 5 HM/2 background for all higher z as well. We stop this run at 2 = 4. 



SPH simulations will be labelled with 'A' and the HYDRA 
ones with 'H'. 



3.1 Gas distribution 

3.1.1 Global distribution 

In figure ^ we present the temperature-density distribution 
for the simulations A-22-64-k and H-22-64-k at z = 2. The 
overall distributions look very similar and can be understood 
in terms of the relative efficiencies of cooling and heating 
and the respective time scales involved, as will be described 
below. We first define the cooling time of gas at temperature 
T as 

" \u\ ^ 2fi p(l-F)2(C-H) ■ ^ ' 

where fcs is Boltzmann's constant, /i is the mean molecular 
weight and tuh is the proton mass. C{u) and TC{p, u) are the 
normalised cooling and heating rates of Appendix |b|. The 
Hubble time is tu = {6-kGp{z))~^^^ . On the left hand panel 
of figure ^ we indicate with a solid line the location where 
the absolute value of the cooling time equals the Hubble 
time at 2 = 2. This line splits the diagram into three sep- 
arate regions. At low densities and high temperatures the 
cooling time is longer than the Hubble time, i.e. neither 
heating nor cooling are able to change the gas tempera- 
ture significantly. At high densities and high temperatures, 
bremsstrahlung and line cooling processes always dominate 
leading to cooling times shorter than the Hubble time. Fi- 
nally, at low densities photo-heating is dominant leading to 
heating times shorter than the Hubble time. (For Ti. < C the 
gas is heated and the cooling time in equation ^ is neg- 
ative. In such cL CcLSG WG C&ll — ^cool 

the heating time, and 

|icooi| the 'net' cooling time.) 

Now consider the line indicated hy Ti. — C, i.e. the equi- 
librium track, at high densities. Approaching this line from 
both lower and higher temperatures we start from a region 
where |fcooi| < and then pass onto the track where the 
denominator of equation (^) tends to zero very fast leading 



to an infinite net cooling time, fcooi 00. Consequently, we 
must go through a point where tcooi ~ tu which explains 
why the solid line is basically identical to the equilibrium 
track at sufficiently high densities. Gas lying just below this 
line will be heated very quickly on a time scale <^ tu onto 
the track, and vice versa for gas slightly hotter than the 
equilibrium T. 

We can now understand the distribution of the gas in 
the {p, T) plane as follows. Efficient photo-heating of under 
dense gas forces it to remain at those temperatures where 
tcooi = tu. Where gas falls into DM potential wells, shock 
heating jgenerates the large plume of hot gas evident in the 
figure. Some of this gas may then reach high enough densities 
so that cooling becomes efficient. This gas condenses onto 
the equilibrium track where heating balances cooling. 

From the previous discussion it is clear that gas at both 
low and high densities is forced by photo-heating to remain 
close to the line where heating is dominant and provides a 
heating (net cooling) time equal to the Hubble time (lower 
branch of the solid line in figure |^) . This line then defines 
a minimum temperature at given density, rmin(p), for the 
simulated gas distribution. Not evident from figured is that 
in fact most of the particles lie very close to this minimum 
temperature. We illustrate this in the following way. We la- 
bel all gas with T < 1.5 x Tmin as being 'cool' and the rest 
as being 'hot' . The condition T — 1.5 x Tmin is shown by 
the dashed line in the figure. The distribution of cool and 
hot gas is shown in detail in figure ^ which compares the 
respective mass fractions as a function of density for vari- 
ous redshifts. The figure shows that at all redshifts most of 
the gas is in the cool phase, though its fraction decreases 
with time. In addition we clearly see an increasing amount 
of gas in the cool phase at the highest densities, resulting 
from cooling in shocked halos at density contrasts > 200. 

Figure ^ shows good agreement between APMSPH and 
HYDRA for the distribution of gas within each phase. The 
distribution of cool gas at low densities (the largest mass 
fraction component) is very similar between APMSPH & 
HYDRA , boding well for Lyman-a simulations. This shows 
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Figure 2. Temperature-density distribution for runs 22-64-k at 2: = 2 for APMSPH (left panel) and HYDRA (right panel). The small 
differences in the low density, low temperature gas distribution between APMSPH and HYDRA are a consequence of the low density 
correction applied to the HYDRA simulation, as discussed in the Appendix. The solid line indicates the locus where the cooling time 
equals the Hubble time, t^^ol = % . At high densities, this corresponds to the equilibrium temperature where cooling balances heating. 
The dashed line shows 1.5xTniin(p), where Tmin is defined as that temperature below which heating is dominant and on which the net 
cooling time is equal to the Hubble time (see text for details). At high densities, T^in equals the equilibrium temperature. Only 1 particle 
in 16 is plotted. 



that the low density correction in HYDRA works well. The 
distribution of hot gas is almost identical between the two 
codes as well, but there are small differences in the high 
density cool phase: APMSPH halos have less cool gas at in- 
termediate densities, but more cool gas at the highest den- 
sities, indicating that APMSPH halos tend to be more con- 
centrated at these highest densities. This is mostly due to 
small differences in the gravitational softening employed in 
these runs (see Table |l|). 

3.1.2 Collapsed systems 

We have run a friends-of-friends group finder with linking 
length 0.177 times the average dark matter inter-particle 
distance on the output files of A-22-64-k and H-22-64-k at 
z — 2.33, thereby selecting halos at an over density of ~ 180. 
Only halos with at least ten dark matter particles are con- 
sidered here. The centre of each dark matter halo is defined 
by the position of the most strongly bound dark dark matter 
particle. We then compute the spherically symmetric den- 
sity profile around the centre of the halo and remove all dark 
matter particles further than r2oo from the centre, thus con- 
fining the halo to its virialised part. Here, r2oo is the distance 
at which the density falls below 200 times the mean density. 
In what follows, SPH particles within h + r2oo from the cen- 
tre of a halo are counted as belonging to that halo. 

The mass function of the halos thus selected is illus- 
trated in figure ^ which shows the number of halos per 
comoving volume as a function of virial mass. There is good 
agreement between the two codes when comparing the total 
masses {i.e. gas plus dark matter) within the virial radius 
'"200- For each halo we have also determined the amount of 



hot cool 




2 4 2 

log (p/Pb) 



Figure 3. Mass fractions at given density in the hot and cool 
phases for various redshifts for runs 22-64-k. Right hand pan- 
els refer to cool gas, left hand panels to hot gas (see text), the 
redshifts are indicated in the different panels. APMSPH results 
are indicated by circles, HYDRA results by squares. Density is in 
units of the average baryon density. The percentages in the panels 
refer to the respective total mass fraction for APMSPH , which 
were very similar for HYDRA . 
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Figure 4. Top panel: mass function for tlie simulation 22-64-k 
at z = 2.33 for APMSPH (full line) and HYDRA (dashed line). 
N denotes the number of halos per comoving Mpc^ with given 
mass in solar units. At this resolution, a halo with mass lO'^'^ Mq 
contains ~ 37 dark matter particles. Bottom panel: same as top 
panel but showing the mass function for the cooled gas in each 
halo. Cooled gas masses are divided by the baryon fraction = 
0.05. 



cooled gas, which for a halo with virial temperature larger 
than 2T4 = 2 X 10*K is gas with temperature less than 2T4,. 
As figure ^ shows, the APMSPH halos tend to have slightly 
more cooled gas than the corresponding HYDRA ones. This 
is at least partly a consequence of the smaller gravitational 
softening in the APMSPH simulations. 

We have computed and compared several other quan- 
tities that characterise the structure and dynamical proper- 
ties of the halos, such as their specific angular momentum 
and their shape parameters, as determined from fitting the 
universal dark matter profiles of Navarro, Frenk & White 
(1996). There is very good agreement between the two codes 
on all these quantities. 

We end this section by showing that the properties of 
the gas distribution within halos are very similar as well. 
Figures (^^) compare the density and temperature distri- 
bution in three halos of total mass within the virial radius of 
5 X 10^^ , 5 X 10" and 5 x lO^^M© respectively. The mean and 
scatter in temperature as well as density at any given radius 
is very similar in the two codes. The APMSPH halos reach 
higher central densities because a smaller gravitational soft- 
ening is used. It is gratifying that the properties of a halo 
with as few as ~ 20 particles of each type are reasonably 
similar when comparing APMSPH with HYDRA . 
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Figure 5. Comparison of gas density (top panels) and temper- 
ature (bottom panels) between APMSPH (left panels) and HY- 
DRA (right panels) for a halo of mass 5 x IO^^Mq in the 22-64-k 
simulations at 2: = 2.33. The APMSPH halo has higher density 
in the centre and also in the in-falling satellites. Thick lines are 
the same in the left and right panels and were drawn to guide the 
eye. The vertical dashed line denotes the gravitational softening, 
the SPH smoothing length h is at least half that. 
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Figure 6. Same as figure ri but for a halo of mass 5 X 10^^ Mq. 



3.2 Spectra 

To compare the statistics of simulated spectra with those 
presented by Hernquist et al. (1996), we will use in this sec- 
tion both their initial conditions and their cooling rates (our 
runs A-22-64-k and H-22-64-k). These simulations should, 
in principle, use the spectral evolution of Haardt & Madau 
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Figure 7. Same as figure H but for a halo of mass 5 X 10^" Mq. 



(1996) but with a reduced amplitude (see Table Q), for which 
we give the fits to the ionising spectrum in Appendix Un- 
fortunately, the actual evolution of the background radiation 
is not identical, since there were small changes between the 
fit to J2i{z) in the preprint of Haardt & Madau (1996), 
from which Hernquist et al. took their evolution, and the 
one which appeared in print (D. Weinberg, private commu- 
nication) . 

This Section is organised as follows. We begin by com- 
paring spectra along a given line of sight between APM- 
SPH and HYDRA . Next, we compare global spectral mea- 
surements (mean absorption and one-point distribution of 
the optical depth). Finally, we compare spectral statistics 
(column- and 6-parameter distributions) deduced from auto- 
mated Voigt-profile fitting of lines. Our detailed descrip tion 
of the calculation of simulated spectra is given in Section A4 



3.2.1 Individual spectra 

We have computed absorption spectra through the middle 
of the computational box of the 22-64-k simulations at a 
redshift of z = 2.33. Figures ^-p^ compare the results of 
APMSPH with HYDRA , showing simulated spectra, den- 
sity, temperature and peculiar velocity along a particular 
line of sight. These figures show that maxima in the den- 
sity distribution give rise to increased absorption with line 
widths determined by Doppler broadening and Hubble ex- 
pansion, as also shown by previous authors. For hydrogen 
there are regions with very low absorption making a distinc- 
tion between 'lines' and continuum possible, but this is no 
longer true for helium nor for hydrogen at higher redshifts. 
The agreement between the two codes is impressive, both 
for Hi and Hen as well as for the total density, over most 
of the velocity range. The excellent agreement increases our 
confidence that both codes are working properly and can be 
used to study the Lyman-a forest. 



Figure 8. Absorption spectrum for Hi through the middle of 
the box of the 22-64-k simulation at z = 2.33. From top to bot- 
tom: Hi absorption spectrum. Hi fraction. Hi weighted tempera- 
ture, Hi weighted peculiar velocity. The top graph is versus veloc- 
ity (wavelength), the bottom three graphs are in units of position 
along the line of sight, from the front of the box to the back. Solid 
lines refer to APMSPH simulation, dashed lines to HYDRA . 



3.2.2 Global spectral characteristics 

We have computed simulated spectra along 1024 lines-of- 
sight on a square grid through the box of the 22-64-k sim- 
ulations at various redshifts. From the spectra we have 
computed the effective mean optical depth at this redshift, 
f{z) = — ln(^^'' e'^ /Np) where the sum is over Np pix- 
els of the spectrum. Figure |l^ shows that APMSPH , HY- 
DRA and TREESPH predict a very similar redshift evo- 
lution for f{z). Most of the remaining differences are pre- 
sumably due to the slightly different assumed evolution of 
the background flux, as can be seen by comparing the HY- 
DRA*and TREESPH* curves. The latter runs use identical 
fits to the photo-ionising fiux although they still use a differ- 
ent fit to the photo-heating rates. A more detailed compar- 
ison is shown in figure |l^ which depicts the one-point prob- 
ability distribution for the optical depth in a pixel. Again 
the comparison is satisfactory, with most of the difference 
between the Tree and the P3M codes being due to the dif- 
ferent background ffux. Note that the apparent shift in the 
mean of P{t) between Tree and P3M codes is not just due to 
the difference in r. Indeed, when the distribution of P{r/f) 
is plotted as in figure |l^ then some differences remain, in- 
dicating that also the shape of P(r) is not quite the same 
between the codes. The APMSPH & HYDRA codes show 
small differences in the high r tail, which is a consequence 
of the small differences in the properties of the high density 
gas, as we remarked on earlier. 



© 0000 RAS, MNRAS 000, 000-000 



10 



T. Theuns et al. 




1 

S 0.5 

0.1 
5 



0.5 



' 1 ' 

o 


APMSPH 


1 1 1 1 I I 1 I 1 ^ 


□ 


HYDRA 




▲ 

- 

' 1 


TREESPH* .^^S^a 


^ \ 


1 o 


1 1 1 1 1 1 1 1 1 1 
APMSPH 


1 [ 1 1 1 1 1 1 r 


□ 


HYDRA 




▲ 


TREESPH* ^..-^^ 










" , 1 







500 1000 

V (km/s) 



1500 2000 



2.5 



3.5 



Figure 9. Same sight line as in figure ^ but now showing 
Hell results. 




Figure 10. Comparison of the state of the IGM between APM- 
SPH (full lines) and HYDRA (dashed line) along the same line 
of sight as figures H and BI showing gas density, temperature and 
peculiar velocity (from top to bottom, respectively) versus posi- 
tion in the box. Density is in units of the average gas density at 
that redshift. 



Figure 11. Mean effective optical depth r as a function of red- 
shift z for runs A-22-64-k & H-22-64-k and published data from 
TREESPH . The top panel refers to hydrogen, the bottom panel 
to singly ionised helium. The APMSPH & HYDRA curves are 
computed from 32^ sight lines through the box. The TREESPH 
data were taken from figure 3 in Croft et al. (1997). The HY- 
DRA*and TREESPH* use the fit equations. 3-5 from Croft et 
al. (1997) for the evolution of the io niza tion rate divided by two, 
the other two use the fit from table p33, divided by two. 




3.2.3 Voigt-profile analysis 

In the previous sections we have shown that APMSPH and 
HYDRA spectra along a given Une of sight look very simi- 
lar and that global statistical quantities of the optical depth 
distribution are similar as well. Here we want to compare the 
statistical properties of the simulated absorption lines based 



Figure 12. One-point probability distribution of optical depth 
for Hi (open symbols) and Hell (filled symbols) for the same sim- 
ulation as shown in figure |ll| at a redshift of z = 2.33. Symbols: 
APMSPH : circles; HYDRA : squares; TREESPH : triangles. The 
APMSPH and HYDRA curves assume the fit in Table ^ for the 
background radiation; the TREESPH data were taken from fig- 
ure 11 in Croft et al. 1997. 
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Figure 13. Same as figure |l^ but for probability distribution of 
the Hi optical depth in units of the mean. 

on Voigt profile (VP) fitting. VP fitting has consistently 
been used to analyse the complex blended features of the 
Lyman-Q forest, seen in real quasar spectra {e.g. Webb 1987, 
Carswell et al. 1987) and in numerical simulations (Dave et 
al. 1997, Haehnelt & Steinmetz 1998, Zhang et al. 1998). 
The technique was devised initially under the premise that 
the forest lines arose from discrete absorbing clouds inter- 
vening the quasar line of sight, a picture currently challenged 
by the successful reproduction of such observed VP parame- 
ter distributions by the smoothly distributed IGM naturally 
occurring within gas-dynamical simulations. 

As we have seen earlier, many of the 'lines' occurring in 
these simulations have significant distortions from the sim- 
ple Voigt profile due to the fact that the structure giving 
rise to the line is extended in space and thus still takes part 
in the Hubble expansion. VP-fitting of such a system then 
requires a complex blend of several Voigt profiles with a 
range in column density and velocity width. Consequently, 
the parameters characterising these components cannot au- 
tomatically be used to infer physical properties of the ab- 
sorbers. Given the loss of a physical justification for the VP 
approach, the method serves merely as way of characterising 
the undulations of what has been described more accurately 
by authors as a Lyman-a Gunn-Peterson-effect (Gunn & 
Peterson 1965). Its power as a diagnostic tool for doing this 
must be scrutinised, before agreement with the observations 
can be described as a success for the model. However, since 
we want to compare our simulated spectra with observed 
ones, it is important that we analyse them in the same way. 
We therefore use a standard package for VP fitting (Carswell 
et al. 1987) used frequently by observers. We will also show 
that the results deduced from VP-fitting are relatively insen- 
sitive to the detailed implementation of the profile fitting. 
There is however one remaining difference between the simu- 
lated and real spectra, namely our simulated spectra are not 
superposed onto the continuum of a background quasar. It is 



therefore difficult to model the continuum fitting performed 
by observers. This problem is aggravated by the small length 
of our simulated spectra compared to observed ones. We try 
to estimate the uncertainty introduced by this below. 

In this paper, we use an automated version of the VP- 
FIT software (Carswell et al. 1987) written to perform 
minimisation over many VP parameters given observational 
spectra. Prior to using VPFIT, each spectrum is convolved 
with a Gaussian profile with FWHM = 8 km s~^, then re- 
sampled onto pixels of width 3 km s~^ to mimic the in- 
strumental profile and characteristics of the HIRES spec- 
trograph on the Keck telescope, which currently provides 
the most up-to-date results on the properties of the weak 
Lyman-a forest lines. Artificial noise is introduced by adding 
a Gaussian random signal with zero mean, and standard de- 
viation a = 0.02 to every pixel (i.e. a SNR of 50 for pixels 
at the continuum), mimicking the read-out noise dominated 
character of modern observed spectra. 

In our simulated spectra, zero-order absorption occurs 
across each line of sight that would presumably be removed 
during observational analysis by the continuum fitting pro- 
cedure. In analysing real spectra, the continuum is normally 
determined by fitting a low-order polynomial to apparently 
'unabsorbed' regions of the spectrum that are typically 
much longer than our simulated spectra. This difference in 
spectral range means that we cannot mimic the continuum 
fitting procedure used by observers on our simulated spec- 
tra. The observational procedure depends a great deal on 
the quality of data as well as on redshift since finding 'un- 
absorbed' regions of the spectrum becomes more difficult 
at higher redshift, where the Lyman-a forest opacities are 
higher. Our spectra cover a small enough velocity range to 
be fit by a fiat continuum, as chosen by a simple procedure 
described as follows. A low average continuum level is as- 
sumed initially, then all pixels below and not within 1 cr of 
this level are rejected, and a new average fiux level for the 
remaining pixels is computed. The same condition is applied 
for this new level and so-on, until the average flux varies by 
less than 1% (note that the signal to noise ratio adopted is 
50 at the continuum). This final average flux level is adopted 
as the fltted continuum, and the spectrum is renormalised 
accordingly. We test for the effects of varying this level on 
the VP results (see below). 

The final stage in preparing spectra for VPFIT is via 
an automatic procedure to find first guess profiles for each 
line. We make initial guesses using the method of Dave et 
al. (1997) for finding weak lines, t.e. firstly smoothing the 
spectrum before 'growing' a Voigt profile into the deepest 
depression. We repeat this procedure on the spectrum with 
the current best fit Voigt profile subtracted, until the resid- 
ual absorption varies by less than 20%. VPFIT turns out to 
be robust in finding accurate fits given few and even quite 
wrong irutial guesses, so we are confident that our results do 
not depend strongly on the details of the procedure for ob- 
taining first guesses. VPFIT then uses these guesses to find 
the best-fitting profiles, adding more lines and removing ill- 
constrained lines where necessary. In general, we calculated 
VP parameters for a spectrum until an overall reduced 
value for the fit falls below 1.3, unless otherwise indicated. 
This is carried out for typically 300 sight line spectra taken 
on a grid through each simulation. We check below the ef- 
fects of varying the required minimum reduced value. 
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It is not clear a priori that the deduced VP statistics 
are independent of the fitting procedure, given e.g. the non- 
uniqueness of Voigt profiles in fitting a blended line and the 
fact that observers often fit profiles by hand. Nevertheless, 
we will show below that our deduced statistics are in excel- 
lent agreement with the ones obtained by Dave et al. (1997), 
using an independent fitting procedure. 

In figure |l^ we show the column density distribution 
(the number A*' of lines per unit redshift and column den- 
sity d?N/ dz /d log A^'hi, or 'differential distribution function' , 
DDF) for runs 22-64-k which use the initial conditions 
from Hernquist et al. (1996) undertaken using APMSPH 
(solid), and HYDRA (dotted) at redshift 2. Also shown are 
the results of a second similar simulation using APMSPH 
(dashed), but where the initial conditions are drawn from a 
different random realization of the initial density field (run 
A-22-64). All three simulations produce very similar DDFs 
across the range of column densities < 10^^ cm~^. There 
is good agreement between the two codes at 2 = 3 as well, 
as illustrated in figure |l^. It is gratifying to conclude that 
these results should not depend greatly on either our code 
implementations, nor the specific realization chosen, nor the 
VP-fitting procedure as applied to two different sets of spec- 
tra at both redshifts. 

Superposed on the simulated DDFs are observation- 
ally determined points of Petitjean et al. (1993, hereafter 
PWRCL) and Kim et al. (1995, hereafter KHCS) for redshift 
z = 2.0 and 2.3, and Hu et al. (1995, hereafter HKCSR), for 
z = 3. We see that the simulations follow the observed points 
well, though possibly under predicting lines with A^hi > lO'^^ 
cm~'^ at 2 = 2. More importantly however, our simulated 
DDFs are directly comparable to the derived column density 
distributions of figures 2 and 3 of Dave et al. (1997) where 
the same cube simulated using TREESPH was analysed in a 
similar way using the AUTOVP/PROFIT fitting software, 
instead of VPFIT. We find that there is a very close similar- 
ity. The results of Dave et al. at 2 = 2 also under predict the 
points of PWRCL and have a very similar slope to our DDF 
for A^Hi > 10^'^ cm~'^. Further at 2 = 3, the distribution lies 
constantly above the 2 = 2 results for Nm > 10^^'^ cm~^, 
but has a shallower slope at lower column densities. The ex- 
act comparison for A^hi < 10^^'^ cm~^, as shown later, is 
sensitive to the VP-fitting conditions, nevertheless the simi- 
larity is striking enough to suggest that the detailed nature 
of the weak Lyman-a forest column density distributions 
are reproducible given both different simulation codes and 
VP-fitting software. 

Figures ([l^p^ show the derived fe-parameter distribu- 
tions for each simulation described above. We follow Dave 
et al. (1997) and include only lines with A^hi > 10^^ cm"^ 
to eliminate broad, weak features, which are sensitive to 
the VP-fitting process, from biasing the comparison with 
observations. These distributions show a consistency indi- 
cating independence with simulation code and realization. 
The shape of the b parameter distribution is very similar to 
the one obtained by Dave et al. (1997, see their figure 3), 
from their analysis of TREESPH simulations with the same 
initial conditions and resolutions as our A-22-64-k and H-22- 
64-k runs. However there is a clear discrepancy between the 
observed b distribution as determined by HKCSR at 2 = 3 
from observations and the simulation results at that red- 
shift. Our simulations all over-predict lines with 6 < 20 and 
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Figure 14. Column density distributions, d? N / dz/ dlogNni, from 
an analysis using VPFIT of 300 spectra along lines of sight 
through several 22Mpc, 64'^ particle simulations, at 2 = 2. Solid 
and dashed lines show results of A-22-64-k, and H-22-64-k simula- 
tions respectively, run using the same initial conditions as used by 
Hernquist et al (1996) (c.f. figure 2 of Dave et al. 1997). The dot- 
ted line shows results of an APMSPH simulation run with differ- 
ent initial conditions (A-22-64). Observational results of PWRCL 
and KHCS with mean absorption redshifts z = 2.0 and 2.3 re- 
spectively, are also plotted as indicated. 

6 > 50 km s~^, whilst failing to reproduce the observed high 
peak of lines at 6 ~ 30 km s~^. Although the 6 parameter 
distribution is sensitive to the signal to noise ratio and the 
VP set-up, we will show below that the dominant reason 
for the discrepancy with observations is lack of numerical 
resolution. 



4 NUMERICAL EFFECTS 

In the previous section we have shown that there is excellent 
agreement between APMSPH , HYDRA and TREESPH on 
a variety of statistics that characterise the simulated Lyman- 
Q spectra. The different codes, however, were compared at 
identical numerical resolution. In this section, we will in- 
vestigate how sensitive these statistics are to details of the 
simulation and analysis, such as numerical resolution, box 
size, and details of the VP-fitting process. To investigate 
the effects of numerical resolution one would ideally like to 
do a simulation of a given region of space at several resolu- 
tions to asses the degree of convergence. However, since the 
required CPU time to perform a given simulation rapidly 
increases with the number of particles in the run, we have 
opted to keep the number of particles fixed but decrease the 
box size. This complicates the analysis since differences in 
statistics between a small and a large box may be due to 
the lack of long wavelength perturbations and saturation of 
modes in the smaller box, rather than non-convergence of 
the simulation in the larger box. In the same vein as the 
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Figure 15. Same as figure |l4| but at 2 = 3. Observational results 
of HKCSR with mean absorption redshifts 2 = 2.9 are plotted. 



Figure 17. b-parameter distributions corresponding to the col- 
umn densities shown in figure Only lines with Af(Hl ) > 10^^ 
cm~^ are included. The histogram shows the corresponding ob- 
servational results of HKCSR. 
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Figure 16. f)-parameter distributions correspondoing to the col- 
umn densities shown in figure Only lines with A'^(Hl ) > 10^^ 
cm~^ are included. 

previous section we will describe first how global properties 
depend on resolution before concentrating on line statistics. 
We will then investigate the dependence of line statistics on 
the VP-fitting procedure. Finally, we will investigate how 
line statistics depend on the assumed UV background. 



4.1 Numerical resolution 



4.1.1 Global spectral characteristics 

We have computed the effective mean optical depth f (2) at 
various redshifts as well as the one-point probability distri- 
bution P(r) from 1024 lines-of-sight for the A-22-64, A-U- 
64 and A-5-64 simulations. Figure hq shows only small differ- 
ences for the thi optical depth as the resolution is increased. 
Note that the optical depth tends to decrease with increasing 
resolution. The main reason for this is that at higher reso- 
lution, small halos collapse that were not resolved at lower 
resolution. This decreases the density of gas in the low den- 
sity regions which leads to a decreasing optical depth. This 
effect is much more pronounced for helium, for which there 
are still significant differences in mean optical depth between 
A- 11-64 and A-5-64. Note that this is a resolution and not 
a box size effect: the A-5-32 box also shown in the figure 
has the same resolution as the A-11-64 box and gives the 
same optical depths despite its different box size. We have 
therefore run an even higher resolution simulation, A-2.5- 
64, which gives a mean optical depth for Hell at 2 = 4 of 
f = 2.63 which is reasonably close to the value 2.78 of the 
A-5-64 simulation. This suggests that the Hell absorption 
has almost converged in our highest resolution simulations, 
but not in the lower resolution ones. 

Zhang et al. (1997) already pointed out that it is far 
more difficult to obtain convergence for the effective mean 
optical depth for helium than for hydrogen. To explain the 
reason for this we first define the effective mean optical depth 
f (p) at given density p from 

exp(-f(p)) = ^exp(-r(p))/iV(p), (3) 

where the sum is over all N{p) pixels in the simulated spec- 
tra where the real space density in the particular pixel is 
p (see e.g. the discussion of figures 6c & 6d in Croft et 
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al. 1997). For later use we define the normalised volume 
fraction of gas at density p as P{p) — N{p)/Np, where Np is 
the total number of pixels. Now, at higher resolution a larger 
fraction of gas collapses into halos with a modest over den- 
sity which are simply not resolved in lower resolution simu- 
lations. At z = 2 for example, this reduces the fraction of gas 
at densities ~ (0.1 — 1) x ps and correspondingly increases 
that fraction at densities ~ (1 — 10) x pB, as is illustrated 
in figure where we sho w the distribution of cool and hot 
gas (as defined in Section 3.1.1 ) at each resolution. The gas 
that collapses to higher densities has an effective mean op- 
tical depth ^ 1 in hydrogen yet r ^ 1 for helium, as is also 
shown in the figure. Therefore, for hydrogen some of the ab- 
sorption lost because of the decreased fraction of low density 
gas is recuperated from the increased opacity in the higher 
density gas. This is not true for helium, however, since the 
lower density gas is already optically thick and increasing its 
column density does not significantly change its net absorp- 
tion. The same reasoning applies at higher redshifts where 
it is even clearer that the amount of cool gas is resolution 
dependent and has only just converged in our highest resolu- 
tion simulation. Miralda-Escude (1993) previously remarked 
that for a large jump in ionising fiux from the Hi limit to the 
Hell limit, as occurs in the Haardt & Madau (1996) spec- 
tra, lines with an Hi column density of ~ 10^^ cm~^ have 
central Hell optical depths ~ 1. Hence for a simulation to 
have a reliable Hell mean effective optical depth it needs to 



resolve well lines with Hi column density of 



10^ 



We will show later (figure |2l|) that the A-22-64 simulation 
produces far fewer lines with such low column densities than 
the higher resolution runs, which then explains why the low 
resolution simulation has a considerably higher mean effec- 
tive Hell optical depth. 

We have shown in figure |l9| that the amount of cool gas 
is strongly resolution dependent. In particular, run 22-64-k, 
which has identical parameters and resolution to the Hern- 
quist et al. (1996) simulation, underpredicts the amount of 
cool, high density gas by a large fraction (a factor of ~ 5 at 
z = 3). Since lines of sight intersecting high density cool gas 
regions produce damped Lyman-a absorbers, the statistics 
of those damped Lyman-a systems, as deduced from these 
simulations (Katz et al. 1996a), are uncertain because they 
are very sensitive to numerical resolution. 

The infiuence of resolution on optical depth is further 
illustrated in figure ^ which compares the net transmission 
of function of its density for our highest and lowest 

resolution runs. The figure illustrates clearly that at higher 
resolution the transmission in the low density gas increases 
significantly but for hydrogen this is mostly compensated 
by a decrease in transmission in the higher density pixels. 
This compensation does not happen for helium and conse- 
quently the helium optical depth decreases noticeably with 
increased resolution. Note also that, although the hydrogen 
mean optical depth has converged there are still significant 
differences in the distribution of P{p) (1 — exp(— f (p))). 



4.1.2 Line statistics 

We have used VP-fitting with the set-up as described earlier 
to characterise the absorption lines in simulated spectra to 
investigate the extent to which VP distributions depend on 
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Figure 18. Mean effective optical dcptti f as a function of red- 
shift 2 for several box sizes indicated in the figure. The top panel 
refers to hydrogen, the bottom panel to singly ionised helium. 
The 64'' simulations have mass resolutions of 1.4 X 10^, 1.8 X 10^, 
2.2 X and 2.1 X Wf' Mq per SPH particle, respectively. 

numerical resolution. The results are shown in figures pl[ - p^ 
for redshifts 2 and 3. The number of hydrogen lines with 
column density > 10^^ cm"'^ is largely independent of reso- 
lution and box size at both redshifts, as shown by the DDFs 
in figures IT] and H, which is quite encouraging. However be- 
low this column density, the lower resolution simulation A- 
22-64 starts to lose a significant number of lines. The higher 
resolution simulations show excellent agreement with the 
observed DDFs also indicated in the figures for both red- 
shifts, although there may be a hint that the simulations 
underproduce lines at redshift 2. 

This agreement does not extent to the 6-parameter dis- 
tributions however. Simulations at different resolutions pro- 
duce quite different 6-parameter distributions, none of which 
fit particularly well the observed one. Higher resolution sim- 
ulations produce a larger fraction of narrower lines. Note 
that the 6-parameter distributions shown in the figures are 
restricted to lines with A^hi > 10^^ cm~'^ where there was 
good agreement for the DDF between simulations at differ- 
ent resolution. We will show in the next section that these 
^-parameter distributions are in fact also quite sensitive to 
the parameters used in the VP- fit analyses (e.ji. required re- 
duced and continuum placement), complicating the com- 
parison with observational data. 



4.2 Dependence on VP-fitting procedure 

The VP-fitting process is sensitive to the quality of the ob- 
servational data, the continuum level chosen, and the strict- 
ness of the fit demanded, as we will quantify here. In fig- 
ure lis] we show the DDF of the high resolution A-11-64 
simulation at z=2, calculated using varying VP-fitting pa- 
rameters. The solid line shows our previous results obtained 
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Figure 19. 
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Mass fractions at given density in the hot and cool phases for var- 
ious redshifts for runs A-22-64 (full lines), A-11-64 (long dashed) 
and A-5-64 (short dashed ) for various redshifts as indicated in the 
panels (see Section |3.1.l| for definition of cool/hot). At higher res- 
olution a larger fraction of gas collapses to higher densities leading 
to less gas at densities p ~ (0.1 — 0.3) X pg which contributes most 
to the absorption. The symbols in the z = b panels refer to our 
highest resolution A-2.5-64 run and compare well with the lower 
resolution A-5-64 run. The dot-dashed lines in the cool gas pan- 
els at 2 = 3 and 2 = 2 indicate the effective mean optical depth 
f (p) as a function of density measured from the simulations for 
hydrogen and helium (curve giving higher values, right scale). 
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Figure 20. Mean transmission at given density, 1 — exp(— f (p)), 
weighted by the fraction by volume of pixels P(p) at that den- 
sity, plotted versus density, illustrating which densities contribute 
most to the absorption as a function of redshift, indicated in the 
panels. Full curves: A-22-64, dashed curves: A-5-64, for hydrogen 
(thick lines) and helium (thin lines). 



using the default reduced requirement, < 1-3, and 
choosing a continuum fit as described previously. The dot- 
ted lines show the results of relaxing the requirement to 
< 2, thereby simulating observations with lower signal to 
noise. This clearly reduces the number of very weak lines, 
A'hi < 10^^'^ cm^, found, but otherwise leaves the DDF un- 
changed. The long-dashed line shows the results when the 
continuum level found by the procedure outlined above is 
then lowered by just 2%, this being the noise level set for 
these spectra (once again requiring < 1.3). Here we see 
a similar but even less noticeable effect. Finally the short- 
dashed line shows the result of raising the continuum level 
by the same amount. We see now that the column density 
distribution is slightly raised down to A^hi < 10^^'^ cm^. 
Beyond this column density it seems that the distribution 
found is insensitive to the chosen parameters of the VP fit- 
ting process. 

Figure ^ shows the 6-parameter results given the dif- 
ferent VP analyses described above. Encouragingly, the 
changes are slight except for the case where the continuum 
level is chosen higher which decreases significantly the height 
of the observed peak and correspondingly increases the num- 
ber of hues at large values of 6 > 100 km s~^. 



4.3 Comparison with UV background scaling 

The deduced DDFs and 6-parameter distributions depend on 
the amplitude of the assumed ionising background which is 
not well known, and on the baryon fraction which is equally 
uncertain. The shape of the radiation spectrum also de- 
termines the amplitude and slope of temperatur e-den sity 
relation obeyed by the cool IGM (see equation (C21) for 
the case of a power law radiation spectrum) and hence has 
some influence on the neutral fraction of gas at a given den- 
sity. This potentially introduces a large parameter space to 
search through with numerical simulations. Fortunately, the 
resulting optical depth only depends on a particular com- 
bination of these parameters. Indeed, in ionization equilib- 
rium the number of photo ionizations per unit time bal- 
ances the number of recombinations, hence J21P oc qhuP^. 
The recombination coefBcient scales with temperature as 
ami oc T^*^'^, so in the low density phase where p oc T^'^ 
we find r oc f2^~°^''^^/"^2i ^ ^^b''^ / J21 {e.g. Hui & Gnedin 
(1997) and Appendix ^ for details). The exponent oi Qb 
depends weakly on the shape of the ionising spectrum. We 
now show that simulations can be run with one set of values 
oi Qb and J21 and later scaled to good accuracy to another 
set. Figure ^ compares the DDFs for simulations H-11-64 
and H-ll-64-j. The first of these was run with the Haardt 
& Madau (1996) ionising background (with amplitude di- 
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Figure 21. Column density distribution dependence on box-size 
and resolution at z=2. Results are shown for three simulations 
run using APMSPH with the same initial modes, but with box 
size and particle numbers as indicated. Observational results of 
KHCS and PWRCL are also plotted as indicated. 



Figure 23. Same as figure |2l| for z = 3. Observational results are 
taken from KHCS and HKCSR, using independent data. 
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Figure 24. Same as figure for z = 3. The histogram shows 
the observational results of HKCSR. 



Figure 22. Corresponding b parameter distributions to the col- 
umn densities shown in figure |2l| for lines with Ni^ > 13 cm~^, 
at 2: = 2. The histogram shows the corresponding observational 
results of KHCS. Note that the highest bin plotted for the KHCS 
results contains only 18 lines. 



vided by 2, indicated by the lines labelled 'HM/2' ) and 
the second with an imposed power law spectrum of ionising 
photons with constant amplitude J21 = 0.5 and slope a = 1. 
The DDFs of these two runs are significantly different, with 
the H- 11-64 run producing more lines at all column densi- 
ties. However, if we assume in the analysts phase that the 
spectrum is HM/2 for the H-ll-64-j run than we obtain an 
almost identical DDF as for the original H- 11-64 run. This 
kind of scaling also works well for the associated 6-parameter 
distributions, shown in figure Esl Note that the run with 
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Figure 25. Column density distributions derived from the same 
A-11-64 simulation a.t z = 2 but varying the reduced require- 
ment and continuum level chosen for each spectrum. Solid and 
dotted lines show results requiring the overall VP-fitting to each 
spectrum to have a reduced < 1-3 S'lid 2 respectively. Short- 
dashed and long-dashed lines show the effect of respectively rais- 
ing and lowering the chosen continuum fit for each spectrum by 
2%, prior to VP-fitting, for x^ < 1.3. Observational results of 
KHCS and PWRCL are also plotted as indicated. 



Figure 27. Column density distributions derived from H-ll-64-j 
simulations at z=2. The solid and dot-dashed lines show results 
where the assumed UV backgrounds are given by HM/2 (tabu- 
lated UV background of HM with amplitude divided by 2), and a 
fixed power-law ( J21 = 0.5, a = 1), respectively both during run- 
time and for the calculation of spectra. Results are also shown for 
the latter simulation but where an HM/2 UV background is as- 
sumed for the calculation of spectra (dashed line), and where the 
original results have been raised by log(2.36) (dotted line). The 
latter factor 2.36 is the ratio of the flux at the Lyman-a limit 
between the J21 = 0.5 power law spectrum and the HM/2 one. 
Observational results of KHCS and PWRCL are also plotted as 
indicated. 
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Figure 26. fe-parameter distributions corresponding to the col- 
umn densities shown in figure ^ for lines with A^hi > lO'^^ cm~^. 
The histogram shows the corresponding observational results of 
KHCS. Lines broader than 100 km s~^ are arbitrarily set to 100 
km s~^. 



power law ionising sources produces a fo-parameter distribu- 
tion very similar to the original HM/2 run, even without 
post-processing scaling. 



5 COMPARISON WITH OBSERVATIONS 

Given the level of convergence of simulated quantities, as 
discussed in the previous section, we now turn our attention 
to a detailed comparison of simulations with observations at 
redshifts > 2. 

A comparison of the effective optical depths for Hi and 
Hell , computed from the A-5-64 run at z=2, with observa- 
tional data is given in figure ^ For the evolution of thi we 
plot the recent results of Fardal et al. (1998), who used a 
combined line-list from a variety of authors for which only 
those originating from studies carried out using the HIRES 
spectrograph on the Keck were included here. We also show 
the effective optical depths observed by Sargent & Barlow 
(published in Ranch et al. 1997). It is clear that the simu- 
lated thi results fit well, as might be expected from the good 
fit of the Hi DDF results. Ranch et al. (1997) used the good 
agreement between simulated and observed optical depths to 
set a lower limit on the baryon fraction, using a lower limit 
on the ionization flux deduced from the observed quasar lu- 
minosity function. Since we confirm their simulation results 
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Figure 28. 6-parameter distributions corresponding to the col- 
umn densities sliown in figure ^ for lines with Afjji > 10^^ cm~^. 
The histogram shows the corresponding observational results of 
KHCS. 

for standard CDM models, we also confirm their lower limit, 
Qb/i^ ~ 0.017. 

For THeii we have used the recent Hell Gunn-Peterson 
effect detections collected by Fardal et al. (1998) (see refer- 
ences therein). Although the available thmi observations are 
limited, the results of Davidsen et al. (1996) do supply a 
strong constraint at Zaha = 2.4. In this case our result that 
the -THcii value inferred from our simulations does indeed 
change significantly when simulating at higher resolution, 
means that this constraint is not matched by our simula- 
tions which predict a thmi lying about a factor 2 below the 
observed value. Unlike the conclusions of Croft et al. (1997) 
we consequently find it is impossible to match both the 
thi and the thou observational constraints by a applying a 
single renormalisation to the Haardt & Madau (1996) UV 
background spectrum, and instead require a 'softer' spectral 
shape. Since thch scales as thch/thi oc rHi/rnoii oc Jhi/Jhcii, 
where Jhi/Jhch is the ratio of the ionising flux at the re- 
spective Hi and Hen limit frequencies, we can increase the 
helium optical depth by scaling the helium ionization rate 
keeping the hydrogen photo-ionization rate constant. In- 
creasing the ratio Jhi/Jhmi by a factor two at all observable 
redshifts from ~ 7 to ~ 14 then provides a good flt to both 
the Hi and Hell optical depths (flgure p9|). 

This softened UV background may yet prove to be 
consistent with the UV background inferred from quasars 
alone, since recent estimates of the 'average' intrinsic quasar 
spectral index have yielded softer values ~ 1.8 (Zheng et 
al. 1997), as opposed to the value 1.5 originally assumed in 
Haardt & Madau (1996). Nevertheless, since the values of 
the UV background at the ionising edges are strongly af- 
fected by the self-absorption by Lyman-Q forest clouds, we 
must wait for more detailed models taking these effects into 
account before drawing any strong conclusions. The values 



Figure 29. Mean optical depth r as a function of redshift z 
for our highest resolution run A-5-64 (open squares connect by 
solid line). The top panel refers to hydrogen, where we have plot- 
ted the observed data points presented in Ranch et al. (1997) 
and those determined by Fardal et al. (1998, denoted FGS) using 
combined line-lists from authors using high-resolution spectra. 
The bottom panel refers to singly ionised helium. The observa- 
tional data points shown are those of Davidsen et al. (1996, DKZ), 
Jakobsen (1997, BDGJP), Reimers et al. (1998, RKWGRW) and 
Tytler et al. (1998, as plotted in figure 1 of Fardal et al. 1998). 
The dashed line connecting triangles in the lower panel shows the 
eff'ect of increasing the FHi/rHeii break in the Haardt &; Madau 
UV background by a factor 2. 

for Fhi, and Fhch depend somewhat on the effective slope of 
the spectrum at frequencies immediately higher than their 
ionising frequencies, however this dependence is weak, so 
achieving the required factor two change in the ratio of 
photo-ionising coefficients by changing the background spec- 
trum in this way is unlikely. 

Previously we have have shown that our SCDM sim- 
ulations using the Haardt & Madau (1996) spectrum (with 
amplitude divided by 2) reproduce the observed DDFs quite 
well, both at z = 2 (figure pi]) and z = 3 (figure The only 
discrepancy may be that the simulations produce a slightly 
lower DDF than seen at z = 2 by Petitjean et al. (1992). 
The difference is slight and in any case the observations are 
as yet not discriminating enough at this redshift to really 
suggest that this is a problem, as the newer results of Kim 
et al. (1997) a.t z — 2.3 are somewhat discrepant with the 
older data. However it is not clear at present to what ex- 
tent the DDF discriminates different plausible cosmological 
models and consequently it is not yet possible to use the 
relatively good agreement between observed and simulated 
DDFs as an argument in favour of our assumed cosmology. 

The simulated and observed b-parameter distributions 
at 2 = 4 are compared in figure ^ We have analysed the 
simulated spectra using a VPFIT set-up that mimics the 
signal to noise ratio and spectral resolution of the data from 
Williger et al. (1994, simulations shown as A-5-64- Williger) 
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Figure 30. fe-paramctcr distribution for lines with A^hi > 10^^ 
cm~^ at 2 = 4 for A-5-64 and A-2.5-64 compared to data from 
Williger et al (1994) and Lu et al. (1996). These two data sets 
were taken with different resolution and so we analysed the sim- 
ulated spectra accordingly: A-5-64- Williger is analysed assuming 
a FWHM of 12 km s"! and SNR of 20 in the VPFIT procedure, 
A-5-64-LU and A-2.5-64-Lu with FWHM of 8 km s"! and SNR 
of 60. 



and Lu et al. (1996, simulations labelled -Lu). For the latter 
data we compare A-5-64 with the higher resolution A-2.5-64 
run. First note the excellent agreement between the diflerent 
resolution runs, suggesting that the 6-parameter distribution 
in the 5Mpc box has very nearly converged. We noted earlier 
that the Hell optical depth of these two runs is very simi- 
lar as well, which increases our confidence in the reliability 
of results drawn from the 5Mpc box run, as far as numer- 
ical artifacts such as resolution are concerned. Next note 
the importance of the wavelength resolution in the analy- 
sis stage, by comparing the A-5-64 run analysed with two 
different VPFIT set-ups. The lower spectral resolution of 
A-5-64- Williger versus A-5-64-Lu reduces dramatically the 
peak in the 6-distribution at ~ 25 km s^^, correspondingly 
increasing the number of broader lines. This trend is also 
shown by the data. The 6-parameter distribution of A-5-64- 
Williger provides an excellent fit to the data from Williger 
et al. (1994), both at low and high h values. At higher reso- 
lution, when we compare the A-5-64-Lu curve with the Lu 
et al. (1996) data, the agreement is still very good. There is 
however a hint that there are too many narrow lines in the 
simulated spectra. 

This difference between simulated and observed distri- 
butions at small &-parameters is even clearer at lower z (see 
figure ^ for z = 2 and figure ^ for z = 3, comparing data 
with the resolved A-5-64 run): the simulated fe-parameter 
distributions peak at lower values of 6 ~ 20 km s^^ whereas 
the observed ones peak at 6 ~ 30 km s~^. A fair compar- 
ison between observations and simulations is hampered by 
the sensitivity of the 6-parameters to the continuum fitting 
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Figure 31. 6-parameter distribution for lines with A'^hi > 10^^ 
cm~^ for A-5-64, after increasing the temperature of all gas arbi- 
trarily by a factor of two, for redshifts 4, 3 and 2 (top to bottom 
respectively). The VPFIT parameters in the analysis were chosen 
to allow comparison with the data of Lu et al. (1996, z = 3.9), 
Hu et al. (1995, z = 2.9) and Kim et al. (1997, z = 2.3), which 
arc superposed for comparison. 



procedure. Note that the high resolution simulations have 
a box size of only ~ lOA, and this might introduce large 
continuum errors. Also, the observed 6-parameter distribu- 
tions have been shown to vary slightly from quasar to quasar 
(Kim et al. 1997), though the effect is larger for the higher 
column density systems A'hi > 10^^ * cm~^. Despite these 
complicating factors the simplest interpretation of our re- 
sults is that our simulated IGM has temperatures that are 
slightly lower by a factor ~ 2 than those existing in the ac- 
tual absorbing IGM. Indeed, if we increase arbitrarily the 
temperature of the simulated gas by a factor of two (dashed 
lines in figure ^l|) then we find excellent agreement between 
the simulated and observed 6-distributions at all redshifts 2, 
3 and 4, in accord with the findings of Haehnelt & Steinmetz 
(1998). 

What could give rise to such a hotter IGM? As dis- 
cussed in Section 2.1, the temperature of the IGM at red- 
shifts z = 4 ^ 2 depends on the epoch of reionization. In 
the Haardt & Madau scenario. Hi reionizes at z ~ 6 and 
Hell at z ~ 4.7 (dashed line in figure lb), and the reioniza- 
tion is due to the increase in the combined QSO luminosity 
at those high redshifts. However, at present there are no 
known QSOs at z > 5, hence the assumed increase in QSO 
flux depends mostly on the extrapolation of the evolution of 
the QSO luminosity function to higher z. If Hell reionization 
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would be delayed, then non-equilibrium efltects would give a 
substantial increase in IGM temperature, even at redshifts 
of « ~ 2. Consequently, uncertainties in the QSO flux at 
redshifts > 5 alone would appear to give sufficient lever- 
age for the required increase in the fGM temperature, ft is 
not clear whether such changes in reionization history alone 
could conspire to give the correct increase in the mean fGM 
temperature over the whole of the redshift interval 4 — > 2. 
Note that any increase in temperature necessary to provide 
fe-distributions that fit observations, would in effect scale 
both the DDF and thi, rneu results further; increasing T by 
a factor of 2 requires increasing fJs by a factor ~ 2°'^ to 
keep the same level of absorption. 

There exist some further numerical possibilities that 
may result in a simulated IGM that appears slightly cool. 
Firstly, larger waves not included in our highest resolution 
run of box-size 5.5Mpc could dynamically heat the medium 
to a greater extent than found in our simulations. The lack 
of dependence of the fo-distribution to the different box-sizes 
examined here indicates that this is probably not an impor- 
tant effect. Secondly, baryons collapsed in small dark matter 
halos that have formed at high-z could possibly be 'evapo- 
rated' by a large temperature boost at reionization such 
as found in the non- equilibrium models. Our smoother tem- 
perature increase during reionization may then allow more 
low-circular velocity halos to capture gas than appropriate. 
To reiterate then, we find that it is likely that a tempera- 
ture boost within our simulations is likely to be necessary 
in order to produce absorption hues that fit the observa- 
tions at « = 2, 3 and 4. The necessary temperature increase 
may well result from a consistent treatment of the thermal 
history as computed by Haardt & Madau (1996) but includ- 
ing non-ionization equilibrium effects, though our estimates 
show that higher temperatures due to an even later reioniza- 
tion epoch, may eventually be necessary. Another plausible 
candidate for extra energy input into the gas is feedback 
from star formation. Certainly it is likely that in future 
observations of the Doppler parameter distribution of the 
Lyman-a forest could become an excellent tool in providing 
constraints on the thermal history of the universe beyond 
« = 2. 



6 SUMMARY AND CONCLUSIONS 

We have presented a new simulation tool, APMSPff , de- 
signed to study numerically the formation of structures re- 
sponsible for the Lyman-a -forest. This code is very fast and 
treats the low density IGM relatively accurately, allowing in- 
creased resolution at little extra simulation time. The IGM is 
allowed to interact with a time-dependent but uniform back- 
ground of ionising photons assumed to come from quasars, 
using the rates suggested by Haardt & Madau (1996). This 
background heats the low density gas and changes the form 
of the cooling function at higher densities. The distribution 
of the gas in the density-temperature plane can be under- 
stood from the relative importance of cooling and heating 
processes, and from the comparison of the appropriate cool- 
ing time scale with the Hubble time. 

We performed extensive comparisons of the new code 
with the HYDRA code of Couchman et al. (1995), which 
was adapted to study this problem as well. The agreement 



between the two codes is excellent for a wide variety of statis- 
tics. The distribution of gas in the (p, T) plane is very similar 
and various statistics on the distribution of halos agree very 
well. The amount of gas which is able to cool in collapsed 
halos is similar in the two codes, showing that coding details 
are not very important in determining this fraction. We are 
currently analysing several large HYDRA simulations per- 
formed on the T3D computer to understand in more detail 
how resolution affects Lyman-a statistics (the VIRGO con- 
sortium, in preparation). Both APMSPH and HYDRA are 
based on the Lagrangian SPH method, which has high res- 
olution in high density regions. However, since many lines 
form in low density regions where SPH suffers from low res- 
olution, it would still be very valuable to compare in detail 
with some of the Eulerian codes used by other groups. 

We also compared our new code with published rc- 
suhs from TREESPH (Hernquist et al. 1996) for simula- 
tions started from their initial conditions and confirm their 
findings. We have also analysed independently our simu- 
lated spectra from these runs using a different implemen- 
tation of automated Voigt profile fitting (VPFIT, CarswcU 
et al. 1987). The deduced line statistics in terms of column 
density and fe-paxameter distributions agree well with their 
published values, showing that Voigt profile fitting gives re- 
producible results. 

We have then used APMSPH to study the effects of 
lack of numerical resolution on quantities deduced from sim- 
ulated spectra based on Voigt profile fitting. The mean ef- 
fective hydrogen optical depth is converged in our medium 
resolution simulation and so are the derived column den- 
sity distributions (DDFs) . The latter also are in good agree- 
ment with DDFs deduced from observations, for our as- 
sumed background fiux and baryon fraction. However, the 
relative amounts of cool gas are rather different when com- 
paring the A-22-64 with the A-11-64 run, which has eight 
times better mass resolution, and there are still noticeable 
differences with our highest A-5-64 run (which has another 
factor of eight better mass resolution) , due to lack of numer- 
ical resolution. With increasing resolution, we find that the 
optical depth decreases, especially for Hell , and that the 
number of lines with small 6-parameter increases. However, 
from a comparison of the A-5-64 run with an even higher 
resolution simulation, A-2.5-64, we find that the A-5-64 box 
is already very close to convergence and we axe relatively 
confident that we can draw reliable conclusions from this 
simulation. We found that the deduced 6-parameter distri- 
butions arc sensitive to the assumed contirmum level, a prob- 
lem which should also infiuence observations to some extent. 
The DDFs, on the other hand, are not very sensitive to the 
exact VP-fitting procedure. 

Some previously published results on the Hen forest are 
unreliable due to lack of numerical resolution. For example 
a,t z — 4, the mean effective Hell optical depths are 4.54, 
3.52, 2.78 and 2.63 for runs A-22-64, A-11-64, A-5-64 and 
A-2.5-64, respectively. This shows that the required resolu- 
tion to get the mean optical depth correct is very high. We 
interpreted the dependence on resolution as being due to 
the formation of progressively smaller halos being resolved 
with better resolution. Low density gas falls into these ha- 
los and hence the optical depth decreases. The good agree- 
ment between the 5.5Mpc and the 2.5Mpc box increases our 
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confidence that these higher resolution runs have effectively 
converged. 

Turning to a comparison of our highest resolution sim- 
ulation with observations, we come to the following conclu- 
sions. 

• There is excellent agreement between the observed and 
simulated Lyman-a column density distributions a.t z = 2 
and 3, provided we divide the ionising background inten- 
sity advocated by Haardt & Madau (1996) by two, for our 
assumed baryon fraction of Qsh^ = 0.0125. Alternatively, 
for the intensity of the ionising background as computed 
by Haardt & Madau, we require a higher baryon fraction 
(Ranch et al. 1997) 



f7s/i~ 0.017 (from DBFs). 



(4) 



• The simulated 6-parameter distributions peak at lower 
6- values than the observed ones for z = 4, 3 and 2, suggesting 
that the simulated IGM temperature in our simulations is 
too low. We argued that uncertainties in reionization history, 
combined with non-equilibrium effects and feedback from 
star formation, might be sufficient to increase the tempera- 
ture by a factor ~ 2, which would bring the simulated dis- 
tributions into excellent agreement with the observed ones. 
However, this would increase the required Q,b even more, 
since increasing the temperature would decrease the amount 
of absorption, giving the higher Q.b limit 



Q,Bh? ~ 0.028 (from 6-parameter distribution) . 



(5) 



• The Hen optical depth corresponding to our best fit 
Hi optical depth is lower than observed values, suggesting 
that the Haardt & Madau ionization spectrum may be too 
hard. The more recent analysis by Zheng et al. (1997) of 
observed quasar spectra lead to a similar conclusion. Fitting 
both Hi and Hell optical depths requires a spectral break 



Jhcu 



14. 



(6) 



Ov prall we find that the level of agreement between aim 



ulat ions of the Lyman - a forest in a scale - imrariant , CDM 
universe and observations, is still impressive. More detailed 
comparisons between simulations and observations will al- 
low us to study the thermal history of the universe at even 
higher redshifts. 
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APPENDIX A: MATHEMATICAL 
DESCRIPTION 

In the first part of this Appendix we give details of the 
equations describing the growth of structure. The second 
part gives the explicit equations for computing SPH quan- 
tities in the APMSPH implementation as well as a detailed 
description of the way we determine SPH neighbours. The 
third part of this Appendix describes the HYDRA low den- 
sity correction. In the final part of this Appendix we give 
the explicit expressions used to compute simulated spectra. 



Al Physical model 

In the Newtonian approximation valid on the scales under 
consideration, the evolution of structures is governed by the 
following set of Lagrangian equations: 



dt ^ 
dx . 1 „ , 1 Vp 

dt a a'' a'^ p 

du „d p 
dt a p 



(Al) 
(A2) 

(A3) 



Here, p = a p is the comoving gas density, Y is the helium 
abundance by mass, (1 — Y) is the hydrogen abundance, 
X = v/a are comoving coordinates, Vp = ax is the peculiar 
velocity, a{t) — {l + z)~^ is the scale factor, t denotes time, z 
is the redshift and V = d/dx. The pressure is p = (7— l)pu, 
where u is the thermal energy per unit mass and 7 = 5/3 for 
a mono-atomic gas. The density pT entering in the Poisson 
equation ( A4 ) is the sum of the dark matter and gas density; 
G is the gravitational constant and niH is the proton mass. 
The dark matter evolves according to the Euler equation 
(A2) with p = 0. These equations neglect feedback from 



stars and AGN. Defining 



K 
U 



P 2 ,3 

-Vpd X 



pud x 



W = / /hpd-^x 



p{H-C)/a^ d^x 



(A5) 
(A6) 
(A7) 
(A8) 



one can write the Layzer-Irvine cosmic energy equation as 
{e.g. Peebles 1980, §24) 



Al = 



{W +{5-3^)0 + a L/d)da 



{K+{3'y-4)U-aL/d)da 

{a - a^)ai{K^ + U^ + W^) 
0, 



(A9) 



where the index i means at the initial expansion factor a;. 

The functions Ti.{p,u) and C{u) in equation (A3) de- 
scribe the heating of the medium by photo-ionization and 
cooling through collisions and interaction with the CMB, 
respectively. In our simulations we use the evolution of the 
photo- ionising background as computed by Haardt & Madau 
(1996). Detailed expressions for the fits to the tempera- 
ture dependence for all included processes are given in Ap- 
pendix [bI 



A2 APMSPH implementation 

The explicit expressions to compute the SPH quantities of 
particle i are 



Mil 

P{i) 



(AlO) 



^ ~ ^TT E^'-^^) - '■(j')) ■ (^(*) - v(j))dw,, 



(All) 



^Vx(i) = ^ ^(x(») - x(j)) ■ (x(*) - x(j)) dW.j 



(A13) 
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Here, Wij = mW (qij) / h'^j is the normalised SPH kernel and 
dWij = mdW{qij)/dqij/qij/hij its derivative; m is the SPH 
particle mass which is the same for all SPH particles. For 
W we use the M4 spline (Monaghan 1992) given by 



W{q) = ^(4- 



6g^ + 3g^) 



if (7 < 1 



= —(2 ~qf if 1 < g < 2 
= else. 



(A14) 



For q < 2/3 we take {l/q)dW{q)/dq = 1/nq in the calcula- 
tion of accelerations to avoid the occurrence of dense knots 
of SPH particles within a gravitational smoothing length. 
We have defined 



qij 



|x(i)-x(j)i 



;ih{i) + Kj)) 



(A15) 
(A16) 
(A17) 



7 

where the artificial viscosity is the sum of a bulk and a von- 
Neumann component: 



H = -a/ipcVrV + /3pft (VrV) 



(A18) 



Note that the latter is in physical (as opposed to comov- 
ing) coordinates: h = ah is the physical smoothing length, 
c = (7(7 — l)^)^'''^ is the physical sound speed and v = 
(d/a)r + ax is the velocity. As is usual, the resolution length 
h is taken such that on average 32 particles ('neighbours' ) 
are within 2h{i) from particle i. For computational efficiency 
we do not allow h to drop below 1/2 the gravitational spline 
softening. We typically take a — 1, P — 2. SPH quantities 
are computed from the particles' positions and velocities in 
two passes over all the neighbours: one pass to co mpute den- 
sity and velocity divergence (equations AlO- All) and a sec- 
ond pass to compute the terms entering in the computation 
of accele rations and thermal energy derivative (equations 



A12 



Ai;) 



Neighbours (particle j is a neighbour of i if its distance 
to i is smaller than twice the SPH smoothing length of i) 
are found using a linked-list (Hockney & Eastwood 1988, 
p. 274). A square grid is placed over the computational vol- 
ume and the linked-list is used as a book-keeping tool to find 
which set of particles resides in a given cell. Since only near- 
est neighbour cells are used to check for potential particle 
neighbours during the SPH loop, to find all neighbours for 
all particles then requires a cell size A = 2/iinax, with /imax 
the maximum smoothing radius of any particle. The usage 
of such a large cell size A becomes rapidly prohibitively ex- 
pensive in evolved systems that have a large dynamic range 
in density and hence in h. We circumvent this problem by 
using a hierarchy of cell sizes: we first loop over all neigh- 
bouring cell pairs using a cell size A but only compute inter- 
actions between particle pairs if at least one of its members 
has //imax < h{i). Typically, we take / = 0.8. In the next 
pass, /imax of those particles for which all forces have not 
been computed yet is now smaller, and so we can do a new 
loop with a smaller cell size for the linked-list until all par- 
ticle pairs have been processed. Consequently, interactions 
between particles in high density regions are computed effi- 
ciently with a small linking cell size yet all potential neigh- 



bours of particles residing in low density regions are still 
found. For systems with a small dynamic range, this extra 
book-keeping leads to a small increase in CPU-time but it 
leads to huge time savings for more clustered systems. In 
fact, the CPU time per step for the SPH calculations in- 
creases only by a factor of 1.6 from a redshift of 50 to the 
final highly clustered redshift of 2 (for a simulation using 
2 X 64^ particles in an L = 22 Mpc box), whereas the CPU 
time for the gravity calculation increases by a factor 3.1 over 
this range. 



A3 HYDRA low density correction 

We now describe the correction to the SPH kernel we em- 
ployed in HYDRA . The effect on the density attributed to a 
particle when the SPH search length is restricted is directly 
dependent on the form of the SPH kernel, W. In fact HY- 
DRA uses exactly the same S PH fo rm a,lisni and kernel as 
APMSPH , so using equations ( AlO ) & (A14) we may write 
the density attributed to any particle as 



E 



(A19) 



where the particle's 'self-density' contribution, = Wa, 
is written explicitly (n is the number of neighbours found 
within the SPH search length 2hi < 2/iinax). We illus- 
trate the results of using this kernel at low densities where 
hi ~ /imax and n is small in figure Al. Here the upper shaded 



region shows the results of using this kernel to find the den- 
sities of 16'' particles randomly filling a box of varying size, 
chosen to represent the range of given baryonic over den- 
sities shown. Each particle's search length, hi, was set to 
the mean inter-particle separation, where this was less than 
/imax . /imax Itsclf was Set at half the mean inter-particle sep- 
aration at over density unity, as this would be the grav- 
itational cell-size imposed during a standard HYDRA run. 
With decreasing p/ps ~ 10, the number of neighbours found 
begins to drop (inset) due to the search length restriction, 
and the calculated particle densities fall slower than do their 
true densities, towards a minimum floor when n = 0. The 
value for the minimum density is just the 'self- density' con- 



tribution specified by the kernel, p" 



/tt in this case. 



Since this minimum density is an arbitrary consequence 
of the kernel used we are essentially free to reset this min- 
imum self-density term as we like. We chose the simplest 
possible modification, calculating densities according to the 
following equation for all particles with n < twenty two 



pi = 



(A20) 



while for all particles with n > 22, the original kernel 
i.e. equation (A19) is still applied. Using this compensated 



kernel, the self-density level is reduced to a value comparable 
to the lowest densities seen in the simulations of APMSPH 
, but the kernel itself is added to, such that for n — 3 1^ th e 



result would be the same as that given by equation (A19). 
We make no changes to the computation of gradients for 
particles with less than 32 neighbours. 

We can now look at the results of using this com- 
pensated kernel to calculate densities of randomly placed 
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gas at difFerent densities as be fore 
lower shaded region of figure Al 



- this is shown as the 
Though the scatter is 
significantly increased, the mean of densities assigned fol- 
lows well the true average densities of the particles down 
to p/pB ~ 0.1. This success is slightly tempered by the 
inevitable dependence of calculated densities on neighbour 
distribution when n is small. For example if the particles 
are not randomly placed, but distributed evenly on a grid, 
then the results of using the cor npen sated kernel to calculate 
particles (dashed line in figure Al) is rather different. The 
calculated densities are biased lower than their true densities 
in a discontinuous fashion since the number n of neighbours 
found falls discontinuously from n = 32 with h < /imax (oc- 
curring at p/pB ~ 10) to 26, 18 and 6 as the search length 
h = /imax, becomes successively smaller in comparison with 
the spacing of the particles at decreasing densities, n finally 
reaches at p/pB = 0, and thus evenly spaced particles 
set up with densities less than this never see any of th eir 
neighbours, and would then according to equation A20 be 
assigned a density simply at the (compensated) self density 
level, po = W i/32 = 8/327r. 

In figure A2 b we show the simulated p~T distribution 
at low densities for one of our runs (H[-22-64-k, see Table ^ 
for run labelling) at z = 2, pre-density reconstruction. Since 
in the simulations the particles are placed initially on a per- 
turbed grid, many particles do initially have n close to zero, 
and are assigned densities biased low towards this compen- 
sated po level, and this situation evidently persists to z—2, 
as there is a group of particles fixed all at the same lowest 
density. However overall the HYDRA p—T distribution com- 
pares very well at lo w de nsities with the distribution found 
by APMSPH (figure A2a) using its exact scheme. This can 
be understood since the evolution of gas in the low density 
IGM with densities biased low is in fact very similar to the 
evolution one would expect with the correct densities. Shock 
heating processes are not underestimated as these are neg- 
ligible for p/pB ~ 10 anyway, and pressure effects are also 
most important for gas settling into DM halos, where the 
gas becomes better and better resolved, and the number of 
neighbours, and hence the accuracy of the assigned densities 
increases. 

We also see that the particles in both APMSPH and pre- 
reconstruction HYDRA lie on a power- law rela tion log(T) oc 



A2|d discussed in 
which 



Qlog(p) (whose form is shown in figure 
detail in Appendix It is this relation which can be 
made use of to post-adjust the particle temperatures at the 
same time as the densities in the reconstruction step. Defin- 
ing pre(post)-reconstruction densities and temperatures as 
Pi(f)jTnf) we fit the slope a obeyed by pi and Ti, find pf ac- 
curately by the density reconstruction step outlined earlier, 
and then set Tf according to 



log(7'/) = logT, -I- a(logp/ - logp,) 



(A21) 



i.e. proportionally adjusting the temperatures along with 
the densities according to the same slope as was obeyed 
pre-reconstruction. We show the results of reconst ruc ting 
the densities and temperatures in this way in figure 



A2 



It 

is clear that this procedure indeed retrieves the low-density 
distribution so that it compares very well with the distri- 
bution simulated by the exact scheme of APMSPH , with 
only a small induced scatter in temperature evident. (The 




log ip/p, 



Figur e A l. Densities calculated us ing k ernels given by equa- 
tions (|A1E|) (upper shaded region) & ( A20 ) (lower shaded region) 



of randomly distributed gas vs true density, where the SPH neigh- 
bour search length, h < /imax, so that the the number of neigh- 
bours found per particle, n (inset) can drop to small values at low 
densities (see text for details). The dashed lines give the results 
for n, and the density calculated using equation (A20) where gas 
is distributed evenly on a grid. 



distribution of gas below the dashed line was shown earlier 
in figure ^. This distribution also compares well with the 
APMSPH one at low densities.) 



A4 Calculation of spectra 

Given the positions, velocities, densities and temperatures 
of all SPH particles at a given redshift, we compute spectra 
along a given line of sight through the box as follows. We 
divide the sight line into A'' ~ 1000 bins of width A in dis- 
tance X along the sight line. For a bin i at position x{i) we 
compute the density and the density weighted temperature 
and velocity from: 



PxU) 



ipT)xij) 



(p'")xij) 



^^E 



X(i)W,, 



(A22) 



(A23) 



X{i){ax{i) + a{x{i) - x{j))W,, 



(A24) 

where X{i) is the abundance of species X of SPH particle i, 
assuming ionization equilibrium {X — }il , X — Hell ; X — 
1 denotes total gas density, see Appendix 



Al 



for meaning 

of other symbols). Labelling bins according to velocity, from 
zero to aL, a bin at velocity v{k) will suffer absorption from 
the bin at velocity v{j) by an amount e~'^''°' where 



r(fc) 



VxU 



-px{j)aA 
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Figure A2. Density-temperature distributions at 2; = 2 of runs 
22-64-k using a) APMSPH , b) HYDRA without reconstruction, 
and c) HYDRA after reconstruction. Inpanel d) the sohd curve 
shows logTuiinCp) for the gas (see fig ure B l) . and its approximated 
power-law form as given by equation (pill). The dotted line shows 
the equilibrium temperature of the gas (where heating balances 
cooling). The dashed lines in each panel are identical and the 
same as those plotted and defined in figure H. 



X exp 



v{k) - v{j) 



Vx{j) 



VUj) = 2kBTx(j)/mx. 



(A25) 



(A26) 



Here, c is the light speed and Vx is the Doppler width of the 
species with particles mass mx- The Lyman-a cross section 
is ac = (37rcrT/8)^''V^o, where ar = 6.625 x 10"^^ cm^ 
is the Thomson cross section, / — 0.41615 is the oscillator 
strength and Ao is the rest wavelength of the transition. For 
the hydrogen Lyman-a transition, we take Ao = 1215. GA 
(a„ =4.45x10"^* cm^), for Hen , Ao = 304. 8 A (ac = 1.12x 
10~^* cm'^. These spectra can converted from 'velocity' v 
to 'observed' wavelength A using A = Ao (1 + z){l + v/c). 



APPENDIX B: ATOMIC PROCESSES AND 
PHOTO-IONIZATION RATES 

Our simulations include all the physical processes relevant 
for the problem of studying primordial gas dynamics in a 
photo ionised intergalactic medium, as first collected by 
Black (1981), and whose completeness is well addressed else- 
where (see e.g. Katz et al. (1996b) and references therein). 
In this Appendix we briefly detail the functional form of 
the atomic physics coefficients used. Both APMSPH and 
HYDRA use the same set of coefficients, all of which are 
based on those collected by Cen (1992) specifically for use in 
cosmological hydrodynamic simulations. A few adjustments 
were made where improvements in accuracy or economy 
were found after comparison with fits used by Efstathiou 
(1992). Where differences occur these have also checked out 



satisfactorily against those quoted using updated atomic 
data by Abel et al. (1997), as detailed below. 

We calculate the normalised radiative cooling function, 
C (normalised as in equation. (A3) such that the rate of 
loss of thermal energy per unit volume, pdu/dt — n'^C) by 
summing the cooling rates whose units and functional de- 
pendence on temperature are given in Table Bl, such that 



C = ^ Ci (T, pb,J{i^,z) 



(Bl) 



The functions used to fit the c; are those of Cen (1992), 
where we have increased the coUisional ionization cooling 
rates by a factor 2 to offset the reducing effects of the 
(1 + (T/lO^Ji')^''^)-^ factor introduced by Cen to extend 
the validity of existing fits to higher temperatures (this fac- 
tor 2 was analogously applied to our coUisional ionization 
rates FeHi.Hci.Hcn given below). 

In order to compute cooling rates as given in Table Bl 
we need ionic abundances for the different gas species. Nor- 
malising fractional densities to the total hydrogen density, 
and denoting them by their standard species nomenclature 
{e.g. Hi = nm/nu), we may write the equations of ioniza- 
tion evolution as 



dHi 
dt 
dHei 

dt 

dHeiii 
di 



ohii WeHii — Hi (Tt-hi + FeHi rie) 

QHoIineHell — Hel (F^Hoi + FeHclJle) 
Hen (F^HcII + FeHoIinc) — «HcIIlHeIII Ue 



supplemented with the closing conditions 

Hi +H11 =1 
Hei -I- Hen -I- Hem = y 
Hll -I- Hen -I- 2Helll = e. 



(B2) 
(B3) 

(B4) 



(B5) 
(B6) 
(B7) 



where y = Y/ {mHe/mH{l — Y)) denotes the helium abun- 
dance by number; TnH,TaH<i are the hydrogen and helium 
atomic mass, ne = enn is the electron number density, and 
the ionization and recombination rates' units and fun ctional 
dependence on temperature are given in Tables B2, B3 & 



Assuming photo-ionization equilibrium we solve the re- 
sulting set of closed equations iteratively until the fractional 
change in all species densities has dropped below 0.001%. 

For the recombination and coUisional ionization rates' 
dependence on temperature T, we used the functional fits 
listed in Table B2. Once again these are based on those of 
Cen (1992), but in addition to the high temperature factor 
correction mentioned above we introduce a factor 0.75 to 
reconcile the fit quoted for the Hll ion recombination rate, 
with the data points given in Spitzer (1978). This inconsis- 
tency between Cen's fit and Spitzer's data has already been 
noted in Rauch et al. (1997, sec. 3.1) where they introduce 
a similar factor 0.8, and does not appear to be present for 
the Hem ion recombination rate coefficient used. 

The rate of photo-ionization of any species ion i, T^i, 
depends on the flux spectrum of ionising ultraviolet back- 



ground photons, J{v, z) (erg cm s Hz sr 



in the 
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following way: 



J{v,z)cr,{v) 
hf 



dv , 



(B8) 



where ai{v) and Vi are the photo- ionization cross-sections 
and ionising threshold frequencies respectively for each 
species. Similarly, excess energy in electrons ejected through 
photo-ionization provides a gas heating mechanism from the 
UV background, for which the normalised pho to-h eating 
function Ti (again normalised as in equation (A3) such 
that the rate of gain in thermal energy per unit volume, 
pdu/dt = n^Ti.) is given by. 



TL = (Hi e^Hi + Hel e^^a + Hell e^Hcn) /"-h ■ 



(B9) 



Here the photo-heating coefficients used, e~f, are calculated 
via an analogous expression to equation 

J(^, z)ai{v)(hv — hv. 



B8 



hv 



■ dv , 



(BIO) 



where hvi is the ionization energy. 

In this paper we assume two separate models for the UV 
background history of the universe. These are implemented 
by deducing photo-ionization and photo-heating rates us- 
ing the above equations. We describe the actual fits imple- 
mented for each coefficient here. 

In case of the evolving background UV spectrum com- 
puted by Haardt & Madau (1996, with deceleration param- 
eter 50 = 0.5, corresponding QSO z-evolution, and QSO 
spectral index equal to 1.5), we performed these integra- 
tions numerically using (Tiiy) given by Cen (1992). Our own 
fits to these integrated photo-ionization and heating rates as 
a function of redshift, given in Table B3, are accurate at all 



redshifts to within 8% and also agree with fits calculated for 
the same spectra independently (Haardt, private communi- 
cation) to within 2% (beyond z = 9 all the rates were set to 
zero). 

As a second more general case we assume a background 
of ionising photons with a power-law spectrum of the stan- 
dard form: 

J{v) = J21 X 10"^' {—\ " ergs-i cm"^ sr"' Hz-\ (Bll) 

where the background fiux is normalised by parameter J21 
at the Hi L yman limit frequency um- We can integrate equa- 
tions [E 



BIO 



in general if we approximate the ion-photon 
cross-sections, Gi, as having simple power-law dependencies 
on frequency. This is already the case for anei, however for 
the hydrogenic ions Hi and Hell we use 



<Ji{v) = 6.3 X 10 



-18 2 ft 

cm 



(B12) 



where Z denotes the atomic number; /i ~ 1 is a dimen- 
sionless constant. By using /hi.hmi ~ 1, 1.21 for the photo- 
ionization rates and /hi.Hch = 1.12, 1.26 for the heating 
rates, the resulting expressions (given in Table B4) agree 



to 2% for a wide range of spectral index, a, with numerical 
integrations using the exact <Ji{v). 



APPENDIX C: EVALUATING THE IGM 
TEMPERATURE 

In this Appendix we write down the defining equations 
needed to solve for the temperature at given density in the 
high redshift universe observable through Lyman- a absorp- 
tion where the net (heating— cooling) time, or 'heating time', 
thcat is equal to the Hubble time tu, 



thcat (p, T, J (u) ,z) = (2) . 



(CI) 



We particularly show that for low densities, this relation ap- 
proximates very well a power-law whose amplitude and slope 
are well-defined, and dependent only weakly on a few param- 
eters, most notably the (effective) UV background spectral 
index (see also Hui & Gnedin 1997). 

We define the Hubble time as follows for an Einstein 
de-Sitter universe. 



1 



4 



(67rGp(2))i/2 (1 + 2)3/2/1 



(C2) 



The heating time (or net cooling time) defined in equa- 
tion (^) may be written in general as 



T 

2(1 - F) Vpo nsABh^Lil + z)3 

^hcat^ 

J7s/l2Ai3L(l + 2)3 ' 



(C3) 



where we use L — Ti. — C to denote the heating rate (erg 
s~^ cm"^), positive for net heating, and A_b = Pb/pb{z). 
Above and in what follows all fundamental atomic parame- 
ters are collected together into constant factors (each given 
a dash), so that the relative contribution from each respec- 
tive process together with cosmological parameters may be 
preserved through the calculation. It should be noted that 
although these factors are denoted by the quantity they at- 
tribute to plus a dash, they will not in general have the same 
units as this quantity. 

we may 



Using the equations as laid out in Appendix Al 



calculate L for a given density and temperature (and UV 
background) and solve equation (Gl) in a numerical fashion 
straightforwardly. We can however simplify the above equa- 
tion considerably by approximating the normalised cooling 
rates in the low-density regime. 

Effectively in the low density region. As < 10, the heat- 
ing rate, L, is dominated by Hi & Hell photo-heating, with a 
small but non-negligible Compton cooling contribution. Fol- 
lowing the notation of the previous Appendix we may write 
the photo- heating component of L, denoted L^, simply as 



em , ehcii 
Hi h Hen 



(C4) 



The Hi Hen and electron fractions above are likely to be 
extremely highly photo-ionised and so we may write their 
abundances accurately as follows 
amine 



Hi 



Hi 



Hen — y 



^ Hcll 



e = l + 2y, 



(C5) 



where the recombination coefficients can be well approxi- 
mated for temperatures T < 10^ K as (see Table B2), 



ami 



aHciii = aHciii-' 



(C6) 



Further, if we model the ph oto-io nization ffux as a power 
law according to equation (Bll), then the ratios of photo- 
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Table Bl. Cooling rates (ergs cm^ s ^). z = redshift, njj = hydrogen number density and T„ = T/(W^K). 
T is in K. 







CoUisional ionization cooling 




Species 


Cl 





2.54 X 10-21rl/2e-157809.1/T(i ^ 




Hi 


C2 


= 


1.88 X 10-21rl/2e-285335.4/T(i ^ 


-Tg/^)-ieHcl 


Hei 


C3 


= 


9.90 X 10-22rl/2e-631515./T(i + 
Recombination cooling 


r5''^)-ieHell 


Hen 


C4 


= 


8.70 X io-27ri/2r3-o-2(i + ro.7) 


-leHll 


Hii 


C5 


= 


1.55 X 10-26r0-3647eHeii 




Hen 


C6 


= 


3.48 X 10-26rl/2T3-0-2(i + 2<>.7) 
Dielectronic recombination coolinj 


"leHelll 


Hem 


C7 




1.24 X io-l3r-l-5e-4™000./T(i ^ 
CoUisional excitation cooling 


- 0.3e-94000/T)gjjgj, 


Heii 


C8 




7.5 X 10-19e-118348/T(i + 2.1/2)- 


-leHl 


Hi 


C9 




5.54 X 10-17T-0-397e-473638/T(i 
Bremsstrahlung 


+ Tg''^)-leHell 


Hen 


CIO 




1.42 X 10-273^Tl/2e(Hll + Hell 


+ 4Helll ) 


Hii , Hell , Hem 


9/ 




1.1 + 0.34e-((5.5-logio(T))2)/3 








Inverse Compton cooling 






Cll 




5.406 X IQ-^^ (T - 2.7(1 + 2))(1 4 







Table B2. Recombination and coUisional ionization rates in s ^, as a function of temperature T(K). 





Recombination 




CKHlI 


= 6.30 X 10-"T-1''2t-0-2/(i + 7-0 




"Hen 


= 1.50 X io-iOt-O-6353 




"Hem 


= 3.36 X 10-iOT-i/2T3-°-2/(1 + Tl 
Dielectronic recombination 


1 ^ "Hcii 


"Hcii 


= 1.9 X 10-3T-i-5e-4-'^xio'/T(i + 


0.3e-9.4xl0 


CoUisional ionization 




TeHl 


= 1.17X 10-10rl/2e-157809.1/T(i. 


^7,1/2)-! 


TeHel 


= 4.76 X 10-llTl/2e-285335-4/T (l . 




rgHcii 


= 1.14 X 10-"Tl/2e-63i5i5/T + 


^1/2)_1 



Table B3. photo-ionization (F^ in s i) and photo-heating (e-^ in ergs i) rates: Haardt & Madau spectrum. 
(z = redshift) 

Ionization rates 

r = e^l+z^2+z^^3 XX X2 X3 

T^m -31.04 2.795 -0.5589 

r.^Hei -31.08 2.822 -0.5664 

r.^Heii -34.30 1.826 -0.3899 

Photo-heating rates 

g^^l + zx^ + z'^X3 ^2 Xs, 

em -56.62 2.788 -0.5594 

e^Hei -56.06 2.800 -0.5531 

e^Heii -58.67 1.888 -0.3947 



Table B4. pho to-iq nization (F-y in s ^) and photo- heating (e-^ in ergs ^) rates: power-law spectra (Compare 
with equation (Bll) for J2i,a definition). 



T-yHel 

r^Hcii 

em 

^7HeII 



Ionization rates 

1.26 X 10-"J2i(3 + a)-l 

1.48 X 10-" J21 0.553" J21 (^jUs 

3.34 X 10-12 j^^ 0.249° (3 + a)-'i 

Photo-heating rates 

2.91 X 10-22 J21 (2 + a)-l(3 + a)-^ 
5.84 X 10- 22 J21 0.553° (^^^ - - 

2.92 X 10-22 J21 0.249° (2 + a)" 



a-l-3.05 ' 



a-f-2.05 ' Q+3.05-' 

1(3-Hq)-i 
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(C7) 



heating and photo-ionization coefficients given in Tabie B4 
are 

CHl _ hpC CHcll _ hpc 

fl^ ~ AH,(2 + a) ' " AHcn(2 + a) ' 

wliere Am and Ahmi are tfie ionization waveiengths for 
Hi and Hen respectiveiy, and we denote Pianck's constant 
and tlie speed of iigfit as hp, c. Drawing ttiese eiements to- 
gether, the normaiised photo-heating rate becomes 



Qhii ^ Q^Hciii 

Am Ahoii 



hpc{l + 2y)T- 
2 + a 



L'T- 



2 + a 



.(C8) 



As mentioned previousiy their is also a smaii contribution 
from Compton cooiing, Lcc, which may be written (cn from 
Table |b^, r> Tcmb) 



5.406x10-^^(1 -h22/)m,i T{l + z) _ L'^^T{l + z) 



(1 - y)po 



nBh^Ap 



.(C9) 



We are now in a position to solve for the temperature 
at a given density where thcat = iff- Substituting for each 



using equations 02 
Qph^AB 



C3 



and 39 



L'T- 



2 + a 



L',,T{l + z ) \ 4(1 + ^ 



\3/2 



■,(C10) 



which by multiplying through by Q,Bh^ Ab /T and juggling 
yields the general solution for T at given A_b, 



T 
To 



tL.th 2 + a 



^cc4(l + 



5/2 



(CU) 



(C12) 



The second factor in brackets encloses the contribution of 
Compton cooling only. This is found always to be small for 
intermediate redshifts. Putting in values for physical con- 
stants, using a cosmological helium fraction, Y = 0.24, and 
^ = (1 + 4y)/(l + y + e) = 0.588 in the low density limit, we 
get 



^Hl 
^hcat 

r' 



2.06x10^'^ s 



(C13) 

5.41x10"" erg cm^ K"^ (C14) 
1.70X 10"^" erg s"^ cm^ K° '' (C15) 
-7.31x10"^'^ erg s"^ cm'* K"^ , (C16) 

where we have used the following atomic numbers for L'^: 

Am = 911.75A (C17) 
Ahcii = 227.67A (C18) 
K'^'^ s-i (C19) 



ami = 2.51x10"" K'^'^ s"^ 



= 1.34x10"^ K" 



(C20) 



which may be collected together into equation ( ^12 ) above 
to find 



To 



3.92xW^{l + z)^{nBh/{2 + a)) — 



K, 



1 + H 



1+^ 



V 5/2\ TT 



h V 1+9.52 y 



(C21) 



where the main dependencies for To are given in the numera- 
tor, and the denominator has the correction due to Compton 
cooling only. 
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